Skip to content

Instantly share code, notes, and snippets.

@ketch

ketch/setplot.py

Created May 29, 2013 14:38
Show Gist options
  • Save ketch/5670771 to your computer and use it in GitHub Desktop.
Save ketch/5670771 to your computer and use it in GitHub Desktop.
Download ZIP
Code for running the LY stegoton experiment with moving wall boundary condition.
Raw
setplot.py
"""
Set up the plot figures, axes, and items to be done for each frame.
This module is imported by the plotting routines and then the
function setplot is called to set the plot parameters.
"""
#--------------------------
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='Stress', figno=1)
# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
#plotaxes.xlimits = [0.,150.]
#plotaxes.ylimits = [-.2,1.0]
plotaxes.title = 'Stress'
# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.plot_var = stress
plotitem.plotstyle = '-o'
plotitem.color = 'b'
plotitem.show = True # show on plot?
plotitem.kwargs = {'linewidth':2,'markersize':5}
# Figure for q[1]
plotfigure = plotdata.new_plotfigure(name='Velocity', figno=2)
# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = 'auto'
plotaxes.ylimits = [-.5,1.1]
plotaxes.title = 'Velocity'
# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.plot_var = velocity
plotitem.plotstyle = '-'
plotitem.color = 'b'
plotitem.show = True # show on plot?
plotitem.kwargs = {'linewidth':3,'markersize':5}
# Parameters used only when creating html and/or latex hardcopy
# e.g., via visclaw.frametools.printframes:
plotdata.printfigs = True # print figures
plotdata.print_format = 'png' # file format
plotdata.print_framenos = 'all' # list of frames to print
plotdata.print_fignos = 'all' # list of figures to print
plotdata.html = True # create html files of plots?
plotdata.html_homelink = '../README.html'
plotdata.latex = True # create latex file of plots?
plotdata.latex_figsperline = 2 # layout of plots
plotdata.latex_framesperline = 1 # layout of plots
plotdata.latex_makepdf = False # also run pdflatex?
return plotdata
def velocity(current_data):
"""Compute velocity from strain and momentum"""
from stegoton import setaux
aux=setaux(current_data.x,rhoB=4,KB=4)
velocity = current_data.q[1,:]/aux[0,:]
return velocity
def stress(current_data):
"""Compute stress from strain and momentum"""
from stegoton import setaux
from clawpack.riemann.rp_nonlinear_elasticity import sigma
aux=setaux(current_data.x)
epsilon = current_data.q[0,:]
stress = sigma(epsilon,aux[1,:])
return stress
Raw
stegoton.py
#!/usr/bin/env python
# encoding: utf-8
import numpy as np
def qinit(state,ic=2,a2=1.0,xupper=600.):
x = state.grid.x.centers
if ic==1: #Zero ic
state.q[:,:] = 0.
elif ic==2:
# Gaussian
sigma = a2*np.exp(-((x-xupper/2.)/10.)**2.)
state.q[0,:] = np.log(sigma 1.)/state.aux[1,:]
state.q[1,:] = 0.
def setaux(x,rhoB=4,KB=4,rhoA=1,KA=1,alpha=0.5,xlower=0.,xupper=600.,bc=2):
aux = np.empty([3,len(x)],order='F')
xfrac = x-np.floor(x)
#Density:
aux[0,:] = rhoA*(xfrac<alpha) rhoB*(xfrac>=alpha)
#Bulk modulus:
aux[1,:] = KA *(xfrac<alpha) KB *(xfrac>=alpha)
aux[2,:] = 0. # not used
return aux
def b4step(solver,state):
#Reverse velocity at trtime
#Note that trtime should be an output point
if state.t>=state.problem_data['trtime']-1.e-10 and not state.problem_data['trdone']:
#print 'Time reversing'
state.q[1,:]=-state.q[1,:]
state.q=state.q
state.problem_data['trdone']=True
if state.t>state.problem_data['trtime']:
print 'WARNING: trtime is ' str(state.problem_data['trtime']) \
' but velocities reversed at time ' str(state.t)
#Change to periodic BCs after initial pulse
if state.t>5*state.problem_data['tw1'] and solver.bc_lower[0]==0:
solver.bc_lower[0]=2
solver.bc_upper[0]=2
def zero_bc(state,dim,t,qbc,num_ghost):
"""Set everything to zero"""
if dim.on_upper_boundary:
qbc[:,-num_ghost:]=0.
def moving_wall_bc(state,dim,t,qbc,num_ghost):
"""Initial pulse generated at left boundary by prescribed motion"""
if dim.on_lower_boundary:
qbc[0,:num_ghost]=qbc[0,num_ghost]
t=state.t; t1=state.problem_data['t1']; tw1=state.problem_data['tw1']
a1=state.problem_data['a1'];
t0 = (t-t1)/tw1
if abs(t0)<=1.: vwall = -a1*(1. np.cos(t0*np.pi))
else: vwall=0.
for ibc in xrange(num_ghost-1):
qbc[1,num_ghost-ibc-1] = 2*vwall*state.aux[1,ibc] - qbc[1,num_ghost ibc]
def stegoton(use_petsc=0,kernel_language='Fortran',solver_type='classic',iplot=0,htmlplot=0,outdir='./_output'):
"""
Stegoton problem.
Nonlinear elasticity in periodic medium.
See LeVeque & Yong (2003).
$$\\epsilon_t - u_x = 0$$
$$\\rho(x) u_t - \\sigma(\\epsilon,x)_x = 0$$
"""
vary_Z=False
if use_petsc:
import clawpack.petclaw as pyclaw
else:
from clawpack import pyclaw
if solver_type=='sharpclaw':
solver = pyclaw.SharpClawSolver1D()
else:
solver = pyclaw.ClawSolver1D()
solver.kernel_language = kernel_language
from clawpack.riemann import rp_nonlinear_elasticity
if kernel_language=='Python':
solver.rp = rp_nonlinear_elasticity.rp_nonlinear_elasticity_1d
elif kernel_language=='Fortran':
from clawpack import riemann
solver.rp = riemann.rp1_nonlinear_elasticity_fwave
solver.bc_lower[0] = pyclaw.BC.custom
solver.bc_upper[0] = pyclaw.BC.extrap
#Use the same BCs for the aux array
solver.aux_bc_lower[0] = pyclaw.BC.extrap#solver.bc_lower
solver.aux_bc_upper = solver.bc_upper
xlower=0.0; xupper=600.0
cellsperlayer=12; mx=int(round(xupper-xlower))*cellsperlayer
x = pyclaw.Dimension('x',xlower,xupper,mx)
domain = pyclaw.Domain(x)
num_eqn = 2
state = pyclaw.State(domain,num_eqn)
#Set global parameters
alpha = 0.5
KA = 1.0
KB = 4.0
rhoA = 1.0
rhoB = 4.0
state.problem_data = {}
state.problem_data['t1'] = 10.0
state.problem_data['tw1'] = 10.0
state.problem_data['a1'] = 0.4
state.problem_data['alpha'] = alpha
state.problem_data['KA'] = KA
state.problem_data['KB'] = KB
state.problem_data['rhoA'] = rhoA
state.problem_data['rhoB'] = rhoB
state.problem_data['trtime'] = 2500.0
state.problem_data['trdone'] = False
#Initialize q and aux
xc=state.grid.x.centers
state.aux=setaux(xc,rhoB,KB,rhoA,KA,alpha,xlower=xlower,xupper=xupper)
qinit(state,ic=2,a2=0.0,xupper=xupper)
tfinal=500.; num_output_times = 10;
solver.max_steps = 5000000
solver.fwave = True
solver.before_step = b4step
solver.user_bc_lower=moving_wall_bc
solver.user_bc_upper=zero_bc
solver.num_waves=2
if solver_type=='sharpclaw':
solver.lim_type = 2
solver.char_decomp=0
claw = pyclaw.Controller()
claw.keep_copy = False
claw.output_style = 1
claw.num_output_times = num_output_times
claw.tfinal = tfinal
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
# Solve
status = claw.run()
if htmlplot: pyclaw.plot.html_plot(outdir=outdir)
if iplot: pyclaw.plot.interactive_plot(outdir=outdir)
if __name__=="__main__":
from clawpack.pyclaw.util import run_app_from_main
output = run_app_from_main(stegoton)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment