Phasefield_students/ 000755 000765 000024 00000000000 13627635253 015512 5 ustar 00benoit staff 000000 000000 Phasefield_students/._.DS_Store 000644 000765 000024 00000000170 13570704534 017404 0 ustar 00benoit staff 000000 000000 Mac OS X 2 F x @ ATTR x x Phasefield_students/.DS_Store 000644 000765 000024 00000014004 13570704534 017170 0 ustar 00benoit staff 000000 000000 Bud1 n b _ c h e . i p y n b _ c h e c k p o i n t slg1Scomp 5š . i p y n b _ c h e c k p o i n t smoDDdutc Ø,›=| . i p y n b _ c h e c k p o i n t smodDdutc Ø,›=| . i p y n b _ c h e c k p o i n t sph1Scomp @ @ € @ € @ € @ E DSDB ` € @ € @ € @ Phasefield_students/._Exercise.html 000644 000765 000024 00000000260 13627635004 020354 0 ustar 00benoit staff 000000 000000 Mac OS X 2 ~ ° ATTR ° com.apple.lastuseddate#PS ':_^ \ÛØ Phasefield_students/Exercise.html 000644 000765 000024 00000233103 13627635004 020143 0 ustar 00benoit staff 000000 000000
Most of the demos of the different equations have been given during the lecture.
We use the following free energy functional: $$ \begin{equation} F=\int f(\phi)+\frac{\alpha}{2}|\nabla \phi|^2 \ dV \label{Free_energy_functional} \end{equation} $$
where the homogeneous free energy density is the simplest double well potential $$ \begin{equation} f(\phi)= 16 \, \Delta f \, \phi^2 \ (1-\phi)^2 \label{dble_well} \end{equation} $$
and where \( \phi \) can be
The double potential is called g
in the python script at the end of the
notebook.
The equilibrium profile of a 1D interface is $$ \begin{equation} \phi_{eq}(x)= \frac12 \left(1+\tanh\left(\frac{2x}{\delta}\right)\right) \label{tanh} \end{equation} $$
where \( \delta=\sqrt{\alpha/(2\Delta f)} \) is the interface width.
The associated interfacial energy is $$ \begin{equation} \gamma = \frac23 \sqrt{2 \,\alpha \, \Delta f} \label{_auto1} \end{equation} $$
Introducing the grid spacing \( d \), the discretized nondimensional free energy \( \tilde{F}={F}/(d^3 \, \Delta f) \) becomes: $$ \begin{equation} \tilde{F} =\sum_n\, \tilde{f}(\phi_n) +\frac{\tilde{\alpha}}{2}|\tilde{\nabla} \phi_n|^2 \label{free_energy_discrete_nondim} \end{equation} $$
where \( \phi_n \) denotes the value of field \( \phi \) at node \( n \), \( n\in\mathbb{N} \), \( \tilde{f}(\phi)= f(\phi)/\Delta f \), \( \tilde{\alpha}=\alpha/(\Delta f\,d^2) \) and \( \tilde{\nabla} \phi \) is the discrete nondimensional gradient operator, for example in 1D \( \phi((n+1)d)-\phi(nd) \).
The nondimensional interface free energy and interfacial width are related to the nondimensional gradient term in the following way: $$ \begin{align} \tilde{\gamma} = \frac{\gamma}{d\,\Delta f} &=\frac{ 2\sqrt{2 }}3 \, \sqrt{\tilde{\alpha}} \label{_auto2}\\ \tilde{\delta} = \frac{\delta}{d} &= \sqrt{\frac{\tilde{\alpha}}{2}} \label{delta_tilde} \end{align} $$
\( \tilde{\delta} \) should be equal to 3 to 5 for numerical reasons (grid pinning).
Assuming a classical relaxation dynamics for the order parameter, the evolution equation reads: $$ \begin{equation} \frac{\partial\phi}{\partial t} = - L \, \frac{\delta F}{\delta \phi} \label{_auto3} \end{equation} $$
where the functional derivative is: $$ \begin{equation} \frac{\delta F}{\delta \phi}= f'(\phi)-\alpha \Delta \phi \label{_auto4} \end{equation} $$
We assume that the characteristic decay time \( t_{exp} \) for a long wavelength fluctuation around \( \phi=0 \) is experimentally known. Using \( f(\phi)\approx \frac{1}{2} f''(0)\,\phi^2 \), and neglecting the heterogeneous term, we have: $$ \begin{equation} \frac{d\phi}{dt} \approx - L \, f''(0) \, \phi \label{_auto5} \end{equation} $$
and finally \( L= 1/ (f''(0) \, t_{exp}) \).
Using the dimensionless time \( \tilde{t}=t/t_{exp} \), the nondimensional Allen-Cahn equation writes: $$ \begin{equation} \frac{d\phi}{d\tilde{t}} = - \tilde L \,\big(\tilde{f}'(\phi)-\tilde{\alpha} \tilde{\Delta} \phi\big) \label{allencahn_nondim} \end{equation} $$
where we have used the kinetic parameter \( \tilde L = {1}/{\tilde{f}''(0)} \) and the discrete laplacian \( \tilde{\Delta}= d^2\Delta \).
Life is often much easier in Fourier space: indeed, differential operators become simple multiplication in reciprocal space. Then, the Allen-Cahn equation in Fourier space writes: $$ \begin{equation} \frac{d\hat\phi}{d\tilde{t}} = - \tilde L \,\left(\,\left\{\tilde{f}'(\phi(\tilde{r}))\right\}_{\tilde{k}} + \tilde{k}^2 \, \tilde{\alpha} \, \hat{\phi} \, \right) \label{eulerspectral} \end{equation} $$
where the Fourier transform of any field \( \psi \) is denoted either as \( \hat{\psi} \) or \( \left\{\psi\right\}_{\tilde{k}} \).
The following explicit first order Euler spectral scheme is already coded: $$ \begin{equation} \hat\phi(\tilde{t}+\tilde{dt}) = \hat\phi(\tilde{t}) - \tilde L \, \tilde{dt} \,\left(\, \left\{\tilde{f}'(\phi(\tilde{r})) \right\}_{\tilde{k},\tilde{t}} + \tilde{\alpha} \,\tilde{k}^2 \,\hat{\phi}(\tilde{t}) \, \right) \label{_auto6} \end{equation} $$
where all quantities on the right hand side are calculated at the previous time step \( \tilde{t} \).
To answer the following questions, several pieces of code are provided
at the end of the present notebook.
a) How can we ensure to have about 3 grid points inside the interface?
Use Eq. \eqref{delta_tilde}. You can check your result by playing with parameter \( \tilde{\alpha} \) considering a plane interface (see next question).
Set alpha
to the right value, and run a calculation considering a planar
precipitate in the middle of the box by commenting out init_interface
with
the proper parameters (and of course commenting all the other functions
changing phi0
).
You should realize that the only driving force is the decrease of the interface energy.
Write down the evolution equation \eqref{allencahn_nondim} in a circular frame, assuming that locally (i.e. in the moving frame of the interface) the profile along the radius remains the same as at equilibrium (this assumption can be checked numerically a posteriori).
Proceed by trial and error: indeed, the formal Von Neumann stability analysis cannot be applied because of the nonlinear term. It may be easier to consider first the case of the lamellar precipitate rather than the circular one. Once a critical value for the time step has been obtained form the lamellar precipitate, check it with the circular one.
The semi-implicit scheme is obtained by evaluating terms linear in \( \hat{\phi} \) on the right hand side of Eq. \eqref{eulerspectral} at time \( t+dt \).
For the next three questions, some calculations must be run on total
time (tmax
) long enough before concluding definitively.
0.01
and 0.1
.
Comment the resulting evolution.
Comment all initial conditions except the two following lines
phi0 = zeros(shape(x_))
phi0 += 0.01*white(shape(x_))
either keeping 0.01
as an amplitude or replacing it by 0.1
.
Comment
phi0 += 0.01*white(shape(x_))
and add at the end of the initial condition (before #20
) the following line
phi0 += 0.5
Don't forget to remove (or comment) the previous add, and comment out the following line:
phi0 = 0.5 + 0.01 * white(shape(x_))
Beware that large driving forces will deviate significantly \( \phi(x) \) across the interface from its equilibrium profile and may give strange and unreliable results.
Determine the value of the driving force necessary to prevent the
dissolution of a circular precipitate of radius \( \tilde{R}^0 = 25 \)
(beware that x.max()
is half the side length of the box).
A planar interface does not move because its energy is invariant with
respect to translations. What happens if the potential \( f \) becomes
non symetric?
Propose a simple way to make it non symetric so as to promote \( \phi=1 \)?
Now \( \phi \) stands for the concentration of some species in a binary alloy.
Consequently, \( \phi \) is conserved (its average value remains constant) and
it must fulfill the following balance law:
$$
\begin{equation}
\frac{\partial\phi}{\partial t} = \nabla\cdot\left(L\nabla\mu\right)
\label{_auto8}
\end{equation}
$$
where \( L \) is the Onsager mobility of the chemical species related to the interdiffusion coefficient, and where \( \mu \) is a generalized chemical potential that reads: $$ \begin{equation} \mu = \frac{df}{d\phi} - \alpha\Delta\phi \label{_auto9} \end{equation} $$
If one assumes that \( L \) does not depend on \( \phi \) (crude approximation) the resulting evolution equation becomes: $$ \begin{equation} \frac{\partial\phi}{\partial t} = L\left(\Delta\frac{df}{d\phi} - \alpha\Delta^2\phi\right) \label{eq:cahnhilliard} \end{equation} $$
Using the same definition for space, time and energy scales as for the Allen-Cahn model, Eq. \eqref{eq:cahnhilliard} can be recast into some nondimensional form (with tildes).
The Cahn-Hilliard equation in Fourier space writes: $$ \begin{equation} \frac{d\hat\phi}{d\tilde{t}} = -\tilde L \,\left(\tilde{k}^2\left\{\tilde{f}'(\phi(\tilde{r}))\right\}_{\tilde{k}} + \tilde{\alpha} \,\tilde{k}^4 \,\hat{\phi} \, \right) \label{eq:cahnhilliard_fourier} \end{equation} $$
It is worth stressing that solving Eq. \eqref{eq:cahnhilliard_fourier} is much easier in Fourier space than in direct space: indeed, the PDE is fourth order in space such that it would require a large discretization stencil to handle in direct space, whereas it is still a simple multiplication in Fourier space.
Using a semi-implicit scheme, we obtain $$ \begin{equation} \hat\phi(\tilde{t}+\tilde{dt}) = \frac{\hat\phi(\tilde{t}) - \tilde L \, \tilde{dt}\, \tilde{k}^2\big\{\tilde{f}'(\phi(\tilde{r})) \big\}_{\tilde{k},\tilde{t}}} { 1 + \tilde L \,\tilde{dt}\,\tilde{\alpha}\,\tilde{k}^4} \label{eq_dphidt_cahnhilliard} \end{equation} $$
a) Code the Cahn-Hilliard equation in a semi-implicit scheme.
You just have to translate in python Eq. \eqref{eq_dphidt_cahnhilliard}.
As for Allen-Cahn, comment all initial conditions except the two following lines:
# Homogeneous at 0.5 with tiny noise
phi0 = 0.5 + 1e-3*white(shape(x_))
Comment the necessary lines in the #(3) Initial microstucture
block.
Comment out any init_circle
and add relevant values to settle
centerx=(0.),centery=(0.)
;radius=(12.)
for an interface width of 5 or 6 pixels
(you can try smaller radii to see what happens);val1=1.
;val2=0.1
for example; you can play with different
supersaturation provided that it remains moderate.
phi0 += init_circle(x_,y_,centerx=(0.),centery=(0.),\
radius=(12.), val1=1., val2=0.1)
'''
Solve Allen-Cahn or Cahn-Hilliard equations using a spectral method
# (The MIT License)
# Copyright (c) 2013-2019
# Benoit Appolaire <benoit.appolaire@univ-lorraine.fr>
# & Yann Le Bouar <yann.lebouar@cnrs.fr>
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
'''
#1
#______________________________________________________________________________
# Import modules/libraries
#______________________________________________________________________________
import sys #system basic routines (may be usefull for debugging)
from scipy import * #scientific library including fast fourier transform
import pylab #plotting a la matlab with matplotlib
import numpy.random as rd #for random numbers
import operator #export functions for the intrinsic operators
#2
#______________________________________________________________________________
# Define a few functions/subroutines
#______________________________________________________________________________
#3
#..............................................................................
def Egrad(phi_,kx_,ky_,alpha):
'''
Gradient energy density
'''
#..............................................................................
tmp = 0.5*alpha*(kx_**2+ky_**2)*abs(phi_)**2
return tmp.real
#4
#..............................................................................
def g(phi):
'''
Double well potential where the energy barrier is 1
'''
#..............................................................................
return 16.*phi**2.*(1.-phi)**2.
#5
#..............................................................................
def dg(phi):
'''
Derivative of the double well potential $g(\phi) = 16*\phi^2(1-\phi)^2$
i.e. $\frac{dg}{dphi}$ (Non linear source term in the pde)
'''
#..............................................................................
return 2.*phi*(1.-phi)*(1.-2.*phi)*16.
#6
#..............................................................................
def ddg(phi):
'''
Second Derivative of the double well potential
$g(\phi) = 16\phi^2(1-\phi)^2$ i.e. $\frac{d^2g}{dphi^2}$
'''
#..............................................................................
return 32.*(1.-6.*phi*(1.-phi))
#7
#..............................................................................
def discretization(L=1.,N=10):
'''
Discretize with N nodes the real and reciprocal spaces
along one direction with length 2*L
'''
#..............................................................................
return mgrid[-L:L:2*L/N], pylab.fftfreq(N,d=L/N)
#8
#..............................................................................
def init_cosx(x_,y_,scalex=1.):
'''
Initialization of a field with cosine along x
'''
#..............................................................................
return (1.+cos(2.*pi*x_/scalex))/2.
#9
#..............................................................................
def init_circle(x_,y_,centerx=0.,centery=0.,radius=1.,val1=1.,val2=0.):
'''
Initialization of a field with a disk centered at (centerx,centery)
and with a given radius
'''
#..............................................................................
return where( (x_-centerx)**2.+(y_-centery)**2.-radius**2.<=0., val1, val2 )
#10
#..............................................................................
def init_interface(x_,y_,pos1=0,pos2=1,val1=1.,val2=0.):
'''
Initialization of a field with two interfaces normal to x
'''
#..............................................................................
return where( (( x_ >= pos1) & ( x_ < pos2)), val1, val2 )
#11
#..............................................................................
def iterwhite():
'''
Generate a sequence of samples of white noise.
Generates a never-ending sequence of floating-point values.
# Copyright (c) 2011 Leif Johnson <leif@leifjohnson.net>
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
'''
#..............................................................................
while True:
for n in rd.randn(100):
yield n
#12
#..............................................................................
def _asarray(source, shape):
'''
# Copyright (c) 2011 Leif Johnson <leif@leifjohnson.net>
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
'''
#..............................................................................
noise = source()
if shape is None:
return noise.next()
count = reduce(operator.mul, shape)
return asarray([noise.next() for _ in range(count)]).reshape(shape)
#13
#..............................................................................
def white(shape=None):
'''
Generate white noise.
shape: If given, returns a numpy array of white noise with this shape.
If not given, return just one sample of white noise.
# Copyright (c) 2011 Leif Johnson <leif@leifjohnson.net>
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
'''
#..............................................................................
return _asarray(iterwhite, shape)
#14
#..............................................................................
def headline():
#..............................................................................
string = '# time'
string += ' fhom grad Ftot '
string += ' dphi <phi> '
string += ' stored in'
print string
#..............................................................................
def monitoring(t,iplot,nl,nc,step,isteps,kx_,ky_,phi0,phi,phi_,alpha):
'''
Print the evolution of global quantities (average, min, max)
and plot snaphots of the phase field
'''
#..............................................................................
fhom = g(phi.real)
grad = Egrad(phi_,kx_,ky_,alpha)
E_grad = real(average(grad))
E_homo = nmax*nmax*average(fhom)
E_tot = E_grad + E_homo
fname = 'phi%02d'%iplot
savez(fname,phi=phi.real)
pylab.subplot(nl,nc,iplot+1)
pylab.imshow(phi.real,interpolation='nearest',origin='lower',\
cmap=pylab.cm.gray,vmin=0.,vmax=1.)
#pylab.contour(phi.real,levels=[0.5],colors='r')
pylab.title('t = %.2f'%(step*pdt))
pylab.xticks([])
pylab.yticks([])
isteps.append(step)
string = '%.2e'%t
string += ' %.3e %.3e %.3e '%(E_grad, E_homo, E_tot)
string += ' %.3e %.3e '%((phi-phi0).real.max(),average(phi.real))
string += ' %s'%(fname+'.npz')
print string
iplot += 1
return iplot
#..............................................................................
def plot_avgminmax(steps,isteps,pdt,tmax,phi_av,phi_min,phi_max):
'''
Plot the evolution of the average, min and max of the phase field.
Can be adapted to plot any other global quantities (e.g. variance)
provided that it has already been computed and pass as an argument.
'''
#..............................................................................
fig2 = pylab.figure(2)
pylab.clf()
pylab.axes([0.14,0.12,0.82,0.83])
trange = array([0]+steps)*pdt
#phi_av = array(phi_av)*100.
pylab.plot(trange,phi_av, 'k-', label='avg')
pylab.plot(trange,phi_min,'b-', label='min')
pylab.plot(trange,phi_max,'r-', label='max')
pylab.legend( loc='center right' )
pylab.axis([0.,tmax,-0.05,1.05])
pylab.xlabel('time', fontsize=20)
pylab.ylabel(r'$\phi$',fontsize=20)
# plot dots corresponding to the snapshots
for i in isteps:
pylab.plot(trange[i],phi_av[i], 'ko',markersize=10)
pylab.plot(trange[i],phi_min[i],'bo',markersize=10)
pylab.plot(trange[i],phi_max[i],'ro',markersize=10)
savetxt('time_table.txt',trange[isteps])
return fig2
#15
#______________________________________________________________________________
# This is the main program which will be run when invoking the script
if __name__ == "__main__":
#______________________________________________________________________________
#This is only for automatic integration by doconce of the
# previous if statement into a jupyter notebook cell to run
dummy_for_ipynb = True
#16
#....................................................
#(0) A few data
#Phase field gradient (nondimensional)
alpha = 1.
#Kinetic coefficient (nondimensional)
kin = 1./ddg(0)
print '# Allen-Cahn mobility = %.2e'%kin
#17
#....................................................
#(1) Space discretization (direct and Fourier)
nmax = 100
xmax = 1.*nmax
x,kx = discretization(L=xmax,N=nmax)
y,ky = discretization(L=xmax,N=nmax)
x_,y_ = meshgrid(x,y)
kx_,ky_ = meshgrid(kx,ky)
kx_ = 2*pi*kx_
ky_ = 2*pi*ky_
#18
#....................................................
#(2) Time discretization
tmax = 1. #total time
pdt = 1e-3
nsteps = int(tmax/pdt)
steps = range(1,nsteps+1)
# Times when to plot
nplots = 10
tplots = linspace(0.,tmax,nplots+1)
#tplots = [80,82,84,86,88,90,92,94,96,98,100,tmax]
iplot = 0
#19
#....................................................
#(3) Initial microstructure: comment
phi0 = zeros(shape(x_))
# Add a vertical slab (lamella) extending from pos1 to pos2
#phi0 += init_interface(x_,y_,pos1=-x.max()/2.,\
# pos2=x.max()/2.,val1=1.,val2=0.)
# Add a circle with given center and radius
#phi0 += init_circle(x_,y_,centerx=(0.),centery=(0.),\
# radius=(x.max()/4.), val1=1., val2=0.)
# Add a circle with given center and radius
#phi0 += init_circle(x_,y_,centerx=(x.max()/2.),\
# centery=(x.max()/6.),radius=(x.max()/3.),\
# val1=1., val2=0.)
# Add noise with 0.01 amplitude
#phi0 += 0.01*white(shape(x_))
# Homogeneous at 0.5 with tiny noise
#phi0 = 0.5 + 1e-3*white(shape(x_))
#20
phi = phi0
phi_ = pylab.fft2(phi)
grad = zeros(shape(x_))
fhom = zeros(shape(x_))
phi_av = [average(phi)]
phi_min = [phi.real.min()]
phi_max = [phi.real.max()]
#21
#....................................................
#(4) Prepare the printing and final plots
headline()
nc = 4
nl = (len(tplots)-1)/nc + 1
fig1 = pylab.figure(1,figsize=(14,3*nl))
step = 0
isteps = []
iplot = monitoring(step*pdt,iplot,nl,nc,step,isteps,\
kx_,ky_,phi0,phi,phi_,alpha)
#22
#....................................................
#(5) Time integration of the PDE
for step in steps:
# (5a) Update the previous field
phi0 = phi
# (5b) Solve in Fourier space using FFT
dg_ = pylab.fft2(dg(phi))
phi_ = pylab.fft2(phi)
#--------------------------------------------------------
#Put here the algebric equation to solve in Fourier space
# Allen-Cahn (explicit time integration scheme)
phi_ = phi_ - kin*pdt*(dg_+alpha*(kx_**2.+ky_**2.)*phi_)
#--------------------------------------------------------
# (5c) Go back to real space
phi = pylab.ifft2(phi_)
# (That's all !!!!!)
# (5d) Compute quantities for the time evolution plot
phi_av.append(average(phi.real))
phi_min.append(phi.real.min())
phi_max.append(phi.real.max())
# (5e) At prescribed steps, print averages and plot snapshots
if step == int(tplots[iplot]/pdt):
iplot = monitoring(step*pdt,iplot,nl,nc,step,isteps,\
kx_,ky_,phi0,phi,phi_,alpha)
#23
#....................................................
#(6) Plot the time evolution of the average and extrema
pylab.subplots_adjust(bottom=0.05,left=0.05,right=0.95,top=0.95,wspace=0.4)
fig2 = plot_avgminmax(steps,isteps,pdt,tmax,phi_av,phi_min,phi_max)
print ; print 'Done' ; print
pylab.show()
# Save
fig1.savefig('snapshots.png')
fig2.savefig('phi_extrema.png')
#24
#....................................................
It can be useful to have a look at a particular snapshot that has been saved in
an npz
file. This can be done using PlotSnapshot.py
which plots side by
side the 2D field and the profile along an horizontal line at position given by
iy
(in nodes).
''' Visualization of npz files generated by the main script
# (The MIT License)
#
# Copyright (c) 2013-2019
# Benoit Appolaire <benoit.appolaire@univ-lorraine.fr>
# & Yann Le Bouar <yann.lebouar@crns.fr>
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
'''
#1
import sys
from pylab import *
from numpy import load
#2
#______________________________________________________________________________
# This is the main program which will be run when invoking the script
if __name__ == "__main__":
#______________________________________________________________________________
name = 'phi10'
res = load(name+'.npz')
phi = res['phi']
iy = 50
#3
figure(1,figsize=(11,5))
ax1 = axes([0.05,0.15,0.4,0.8])
imshow(phi,interpolation='nearest',origin='lower',\
cmap=cm.gray,vmin=0.,vmax=1.)
axhline(iy,color='r')
xticks([])
yticks([])
#4
ax2 = axes([0.53,0.15,0.44,0.8])
plot(phi[iy,:],'r-o',linewidth=1.5)
axis([0,len(phi)-1,-0.01,1.01])
xlabel('$x$',fontsize=24)
ylabel('$\phi$',fontsize=24)
savefig(name+'.png')
show()
#5