# 2-dimensional p-system¶

## Two-dimensional p-system¶

Solve the two-dimensional generalization of the p-system:

$\begin{split}\epsilon_t - u_x - v_y & = 0 \\ \rho(x,y) u_t - \sigma(\epsilon,x,y)_x & = 0 \\ \rho(x,y) v_t - \sigma(\epsilon,x,y)_y & = 0.\end{split}$

We take $$\sigma = e^{K(x,y)\epsilon} - 1$$, and the material coefficients $$\rho,K$$ vary in a checkerboard pattern. The resulting dynamics lead to solitary waves, though much more resolution is needed in order to see them.

This example shows how to set an aux array, use a b4step function, use gauges, compute output functionals, and restart a simulation from a checkpoint.

## Source:¶

#!/usr/bin/env python
# encoding: utf-8
r"""
Two-dimensional p-system
==============================

Solve the two-dimensional generalization of the p-system:

.. math::
\epsilon_t - u_x - v_y & = 0 \\
\rho(x,y) u_t - \sigma(\epsilon,x,y)_x & = 0 \\
\rho(x,y) v_t - \sigma(\epsilon,x,y)_y & = 0.

We take :math:\sigma = e^{K(x,y)\epsilon} - 1, and the
material coefficients :math:\rho,K vary in a checkerboard
pattern.  The resulting dynamics lead to solitary waves,
though much more resolution is needed in order to see them.

This example shows how to set an aux array, use a b4step function,
use gauges, compute output functionals, and restart a simulation
from a checkpoint.
"""

from __future__ import absolute_import
import numpy as np
from clawpack import riemann
from six.moves import range

def qinit(state,A,x0,y0,varx,vary):
r""" Set initial conditions:
Gaussian stress, zero velocities."""
yy,xx = state.grid.c_centers
stress=A*np.exp(-(xx-x0)**2/(2*varx)-(yy-y0)**2/(2*vary)) #sigma(@t=0)
stress_rel=state.aux[2,:]
K=state.aux[1,:]

state.q[0,:,:]=np.where(stress_rel==1,1,0)*stress/K+np.where(stress_rel==2,1,0)*np.log(stress+1)/K
state.q[1,:,:]=0; state.q[2,:,:]=0

def setaux(x,y, KA=1, KB=4, rhoA=1, rhoB=4, stress_rel=2):
r"""Return an array containing the values of the material
coefficients.

aux[0,i,j] = rho(x_i, y_j)              (material density)
aux[1,i,j] = K(x_i, y_j)                (bulk modulus)
aux[2,i,j] = stress-strain relation type at (x_i, y_j)
"""
alphax=0.5; deltax=1.
alphay=0.5; deltay=1.

medium_type = 'checkerboard'

aux = np.empty((4,len(x),len(y)), order='F')
if medium_type == 'checkerboard':
# xfrac and yfrac are x and y relative to deltax and deltay resp.
xfrac=x-np.floor(x/deltax)*deltax
yfrac=y-np.floor(y/deltay)*deltay
# create a meshgrid out of xfrac and yfrac
[yf,xf]=np.meshgrid(yfrac,xfrac)
# density
aux[0,:,:]=rhoA*(xf<=alphax*deltax)*(yf<=alphay*deltay)\
+rhoA*(xf >alphax*deltax)*(yf >alphay*deltay)\
+rhoB*(xf >alphax*deltax)*(yf<=alphay*deltay)\
+rhoB*(xf<=alphax*deltax)*(yf >alphay*deltay)
#Young's modulus
aux[1,:,:]=KA*(xf<=alphax*deltax)*(yf<=alphay*deltay)\
+KA*(xf >alphax*deltax)*(yf >alphay*deltay)\
+KB*(xf >alphax*deltax)*(yf<=alphay*deltay)\
+KB*(xf<=alphax*deltax)*(yf >alphay*deltay)
# linearity of material
aux[2,:,:]=stress_rel
elif medium_type == 'sinusoidal' or medium_type == 'smooth_checkerboard':
[yy,xx]=np.meshgrid(y,x)
Amp_rho=np.abs(rhoA-rhoB)/2; offset_p=(rhoA+rhoB)/2
Amp_K=np.abs(KA-KB)/2; offset_E=(KA+KB)/2
if medium_type == 'sinusoidal':
frec_x=2*np.pi/deltax; frec_y=2*np.pi/deltay
fun=np.sin(frec_x*xx)*np.sin(frec_y*yy)
else:
sharpness=10
fun_x=xx*0; fun_y=yy*0
for i in range(0,1+int(np.ceil((x[-1]-x[0])/(deltax*0.5)))):
fun_x=fun_x+(-1)**i*np.tanh(sharpness*(xx-deltax*i*0.5))
for i in range(0,1+int(np.ceil((y[-1]-y[0])/(deltay*0.5)))):
fun_y=fun_y+(-1)**i*np.tanh(sharpness*(yy-deltay*i*0.5))
fun=fun_x*fun_y
aux[0,:,:]=Amp_rho*fun+offset_p
aux[1,:,:]=Amp_K*fun+offset_E
aux[2,:,:]=stress_rel
return aux

def b4step(solver,state):
r"""Put in aux[3,:,:] the value of q[0,:,:] (eps).
This is required in rptpv.f.
Only used by classic (not SharpClaw).
"""
state.aux[3,:,:] = state.q[0,:,:]

def compute_stress(state):
""" Compute stress from strain and store in state.p."""
K=state.aux[1,:,:]
stress_rel=state.aux[2,:,:]
eps=state.q[0,:,:]
state.p[0,:,:] = np.where(stress_rel==1,1,0) * K*eps \
+np.where(stress_rel==2,1,0) * (np.exp(eps*K)-1) \

def total_energy(state):
rho = state.aux[0,:,:]; K = state.aux[1,:,:]

u = state.q[1,:,:]/rho
v = state.q[2,:,:]/rho
kinetic=rho * (u**2 + v**2)/2.

eps = state.q[0,:,:]
sigma = np.exp(K*eps) - 1.
potential = (sigma-np.log(sigma+1.))/K

dx=state.grid.delta[0]; dy=state.grid.delta[1]

state.F[0,:,:] = (potential+kinetic)*dx*dy

def gauge_stress(q,aux):
p = np.exp(q[0]*aux[1])-1
return [p,10*p]

def setup(kernel_language='Fortran',
use_petsc=False,outdir='./_output',solver_type='classic',
disable_output=False, cells_per_layer=30, tfinal=18.):

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

# material parameters
KA=1.;   rhoA=1.
KB=4.;   rhoB=4.
stress_rel=2;

# Domain
x_lower=0.25; x_upper=20.25
y_lower=0.25; y_upper=20.25
# cells per layer
mx=int((x_upper-x_lower)*cells_per_layer)
my=int((y_upper-y_lower)*cells_per_layer)
# Initial condition parameters
initial_amplitude=10.
x0=0.25 # Center of initial perturbation
y0=0.25 # Center of initial perturbation
varx=0.5; vary=0.5 # Width of initial perturbation

num_output_times = 10

if solver_type == 'classic':
solver = pyclaw.ClawSolver2D(riemann.psystem_2D)
solver.dimensional_split=False
solver.cfl_max = 0.9
solver.cfl_desired = 0.8
solver.limiters = pyclaw.limiters.tvd.superbee
elif solver_type=='sharpclaw':
solver = pyclaw.SharpClawSolver2D(riemann.psystem_2D)

if kernel_language != 'Fortran':
raise Exception('Unrecognized value of kernel_language for 2D psystem')

# Boundary conditions
solver.bc_lower[0] = pyclaw.BC.wall
solver.bc_upper[0] = pyclaw.BC.extrap
solver.bc_lower[1] = pyclaw.BC.wall
solver.bc_upper[1] = pyclaw.BC.extrap
solver.aux_bc_lower[0] = pyclaw.BC.wall
solver.aux_bc_upper[0] = pyclaw.BC.extrap
solver.aux_bc_lower[1] = pyclaw.BC.wall
solver.aux_bc_upper[1] = pyclaw.BC.extrap

solver.fwave = True
solver.before_step = b4step

#controller
claw = pyclaw.Controller()
claw.tfinal = tfinal
claw.solver = solver
claw.outdir = outdir

# restart options
restart_from_frame = None

if restart_from_frame is None:
x = pyclaw.Dimension(x_lower,x_upper,mx,name='x')
y = pyclaw.Dimension(y_lower,y_upper,my,name='y')
domain = pyclaw.Domain([x,y])
num_eqn = 3
num_aux = 4
state = pyclaw.State(domain,num_eqn,num_aux)
state.mF = 1
state.mp = 1

grid = state.grid
state.aux = setaux(grid.x.centers,grid.y.centers,KA,KB,rhoA,rhoB,stress_rel)
#Initial condition
qinit(state,initial_amplitude,x0,y0,varx,vary)

claw.solution = pyclaw.Solution(state,domain)
claw.num_output_times = num_output_times

else:
claw.solution = pyclaw.Solution(restart_from_frame, format='petsc',read_aux=False)
claw.solution.state.mp = 1
grid = claw.solution.domain.grid
claw.solution.state.aux = setaux(grid.x.centers,grid.y.centers)
claw.num_output_times = num_output_times - restart_from_frame
claw.start_frame = restart_from_frame

#claw.p_function = p_function
if disable_output:
claw.output_format = None
claw.compute_F = total_energy
claw.compute_p = compute_stress
claw.write_aux_init = False

solver.compute_gauge_values = gauge_stress
state.keep_gauges = True
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.
"""
from clawpack.visclaw import colormaps

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

# Figure for strain
plotfigure = plotdata.new_plotfigure(name='Stress', figno=0)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = 'Strain'
plotaxes.xlimits = [0.,20.]
plotaxes.ylimits = [0.,20.]
plotaxes.scaled = True

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

return plotdata

def stress(current_data):
import numpy as np
aux = setaux(current_data.x[:,0],current_data.y[0,:])
q = current_data.q
return np.exp(aux[1,...]*q[0,...])-1.

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