Phase field: practical session

B. Appolaire

Université de Lorraine, France

March 2020


Most of the demos of the different equations have been given during the lecture.


Free energy functional


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

Warning.

The double potential is called g in the python script at the end of the notebook.




Problem 1: Allen-Cahn model


Equilibrium properties

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).


Evolution equation

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 \).


Solving in Fourier space

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} \).



Questions

Warning.


To answer the following questions, several pieces of code are provided at the end of the present notebook.


Beware to indent using spaces rather than tabs (one level = 3 spaces) when modifying the python scripts.


a) How can we ensure to have about 3 grid points inside the interface?

Hint.

Use Eq. \eqref{delta_tilde}. You can check your result by playing with parameter \( \tilde{\alpha} \) considering a plane interface (see next question).


b) Compute an equilibrium flat profile and check that it is compliant with the theoretical profile. Compute the nondimensional interfacial free energy. Compare it to the theoretical value.

Hint.

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).


c) Start a simulation with a circle and study the evolution of the average phase field. Explain the result.

Hint.

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).


d) What is the maximum time step that can be used?

Hint.

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.


e) Code a semi implicit spectral scheme and analyze its stability and accuracy (you may use the simulation of a shrinking circle).

Hint.

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 \).

Warning.

For the next three questions, some calculations must be run on total time (tmax) long enough before concluding definitively.


f) Start the simulation with \( \phi=0 \) everywhere plus some noise with amplitude 0.01 and 0.1. Comment the resulting evolution.

Hint.

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.


g) Start the simulation with \( \phi=0.5 \) everywhere. Comment the evolution.

Hint.

Comment

   phi0   += 0.01*white(shape(x_))

and add at the end of the initial condition (before #20) the following line

   phi0 += 0.5


h) Start the simulation with \( \phi=0.5 \) everywhere plus some tiny noise. Comment the evolution.

Hint.

Don't forget to remove (or comment) the previous add, and comment out the following line:

   phi0 = 0.5 + 0.01 * white(shape(x_))


i) What should be modified to promote the growth of regions with \( \phi=1 \)? Implement this modification and deduce the growth law for a lamella.

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).

Hint.

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 \)?




Problem 2: Cahn-Hilliard model


Evolution equation

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).


Solving in Fourier space

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} $$


Questions


a) Code the Cahn-Hilliard equation in a semi-implicit scheme.

Hint.

You just have to translate in python Eq. \eqref{eq_dphidt_cahnhilliard}.


b) Run a calculation starting with an initial concentration at \( \phi=0.5 \) plus some tiny noise. Compare qualitatively the long time evolution with that of Allen-Cahn with the same initial conditions.
Play with similar initial conditions but different average values of \( \phi \). Compare the different evolutions.

Hint.

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_))


c) Settle some clever initial condition to investigate the growth of a small circular rich precipitate (the initial radius should be at least twice the interface width).

Hint.

Comment the necessary lines in the #(3) Initial microstucture block. Comment out any init_circle and add relevant values to settle

  • the disk center at the center of the box centerx=(0.),centery=(0.);
  • the radius of 12 pixels radius=(12.) for an interface width of 5 or 6 pixels (you can try smaller radii to see what happens);
  • a concentration of 1 in the precipitate val1=1.;
  • and a concentration out of equilibrium in the matrix (supersaturation here is equal to the concentration of the matrix because \( \Delta c^{\rm{eq}}=1 \) in this simple problem) val2=0.1 for example; you can play with different supersaturation provided that it remains moderate.
Hence, the initial conditions look like this:

   phi0   += init_circle(x_,y_,centerx=(0.),centery=(0.),\ 
                         radius=(12.), val1=1., val2=0.1)


d) Determine the growth law of the circular precipitate. Comment.


e) Settle a simple but clever initial condition to investigate Ostwald ripening process that play the major role in coarsening microstructures when holding a two-phase microstructure at high temperature.



The python scripts


''' 
 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
   #....................................................



Plotting a snapshot

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

© 2019, B. Appolaire. CC BY-NC-SA