Lecture5.tex

\documentclass[pdf,ps,8pt]{beamer}

\usepackage{lmodern}
\usepackage{amsmath}
\usepackage{color}
\usepackage{bbold}
\usepackage{cancel}
\usepackage{slashed}
\usepackage{graphicx}
\usepackage{verbatim}
\usepackage{textcomp}

\usetheme{Singapore}

\definecolor{palegray}{rgb}{0.82,0.822,0.82}
\newcommand{\preliminary}{
{ \rput{30}(7,-4.0){\fontsize{40}{40}\selectfont {\color{palegray}Preliminary Preliminary}} }
}

\newcommand{\textapprox}{\raisebox{0.5ex}{\texttildelow}}

\newcommand{\miniscule}{\fontsize{3pt}{4pt}\selectfont}

\def\MSbar{$\overline{\mathrm{MS}}$}
\def\gev{\,\mathrm{GeV}}
\def\mev{\,\mathrm{MeV}}
\def\fm{\,\mathrm{fm}}
\def\SU{\mathrm{SU}}
\def\su#1#2{\SU(#1)_\mathrm{#2}}
\def\rpisq{\langle r_\pi^2\rangle}
% final values
\def\rpisqsim{0.38(4)}  % pole fit for 330 MeV pion
\def\rpisqsu2sim{0.354(31)}  % su2 fit for 330 MeV pion
\def\rpisqsimlong{0.382(42)}
\def\rpisqphys{0.418(31)} % SU(2) chi extrap
\def\rpisqphyslong{0.418(31)}
\def\lsixr{-0.0093(10)} % SU(2)
\def\Lniner{0.0031(6)}  % SU(3)

\newcommand{\chpt}{\chi^{\rm PT}}
\newcommand{\tchpt}{$\chi^{\rm PT}$}
\newcommand{\tchptthree}{$\chi^{\rm PT}_3$}
\newcommand{\tchpttwo}{$\chi^{\rm PT}_2$}

\newcommand{\xiav}{\langle\,\xi\,\rangle}
\newcommand{\xisqav}{\langle\,\xi^2\,\rangle}


\newcommand{\mD}{\left(\begin{array}{cc} \DO & \Dd \\  \Ddb&\DOb \end{array} \right)}

\newcommand{\Ob}{\bar{\Omega}}

\newcommand{\DO}{D_\Omega}
\newcommand{\Dd}{D_\partial}
\newcommand{\DOi}{D_\Omega^{-1}}
\newcommand{\DOid}{D_\Omega^{-\dagger}}
\newcommand{\Pd} {\mathbb{P}_\partial}
\newcommand{\PO} {\mathbb{P}_\Omega}

\newcommand{\DOb}{D_{\bar{\Omega}}}
\newcommand{\Ddb}{D_{\bar{\partial}}}
\newcommand{\DObi}{D_{\bar{\Omega}}^{-1}}
\newcommand{\DObid}{D_{\bar{\Omega}}^{-\dagger}}
\newcommand{\Pdb}{\mathbb{P}_{\bar{\partial}}}
\newcommand{\POb} {\mathbb{P}_{\bar\Omega}}

\newcommand{\Phidb}{\mathbb{\phi}_{\bar{\partial}}}
\newcommand{\etadb}{\mathbb{\eta}_{\bar{\partial}}}

\newcommand{\hDO}{\hat D_\Omega}
\newcommand{\hDd}{\hat D_\partial}
\newcommand{\hDOi}{\hat D_\Omega^{-1}}
\newcommand{\hPd} {\hat{\mathbb{P}}_\partial}

\newcommand{\hDOb}{\hat D_{\bar{\Omega}}}
\newcommand{\hDdb}{\hat D_{\bar{\partial}}}
\newcommand{\hDObi}{\hat D_{\bar{\Omega}}^{-1}}
\newcommand{\hPdb}{\hat{\mathbb{P}}_{\bar{\partial}}}

\newcommand{\mul}[1]{\left(\begin{array}{cc}#1 & 0 \\ 0& 0\end{array}\right)}
\newcommand{\mur}[1]{\left(\begin{array}{cc}0  & #1\\ 0& 0\end{array}\right)}
\newcommand{\mll}[1]{\left(\begin{array}{cc}0  & 0 \\ #1 & 0\end{array}\right)}
\newcommand{\mlr}[1]{\left(\begin{array}{cc}0  & 0 \\ 0& #1\end{array}\right)}

\newcommand{\mDO}{\mul{ \DO}}
\newcommand{\mDd}{\mur{ \Dd}}
\newcommand{\mDOi}{\mul{\DOi}}
\newcommand{\mPd} {\mlr{\Pd}}

\newcommand{\mDOb}{\mlr{\DOb}}
\newcommand{\mDdb}{\mll{\Ddb}}
\newcommand{\mDObi}{\mlr{\DObi}}
\newcommand{\mPdb}{\mul{\Pdb}}
\newcommand{\rmod}{\mathrm{mod}}
\newcommand{\rdiv}{\mathrm{div}}

\newcommand{\link}[1]{\href{#1}{ {\color{blue} #1} }}

\beamertemplatenavigationsymbolsempty
\begin{document}

\begin{frame}[fragile]\small\frametitle{  Computational Methods (practice) -  Lecture 5    }

  \begin{center}
 
  {\color{red} Peter Boyle} (BNL, Edinburgh)

  \begin{itemize}
  \item Walk through an example code: u/d/s/c Meson spectrum
  \item Under 300 lines of code
  \item How to make Point, Z2, Gaussian, Sequential propagators
  \item Final words
  \end{itemize}

\end{center}  

\end{frame}


\begin{frame}[fragile]\small\frametitle{Observables}

Importance sampling has reduced:
$$
\langle {\cal O} \rangle = \frac{1}{Z} \int_U e^{-S_G[U]} {\cal O}(U) dU  \to \frac{1}{N} \sum\limits_i {\cal O}(U_i) 
$$
\begin{itemize}
\item Zero momentum pion, kaon or B meson two point function
  $$
  \sum\limits_x \langle \bar u \gamma_0 \gamma_5 d (x,t) \bar d \gamma_0 \gamma_5 u (0,0) \rangle
    = \frac{1}{N}\sum_i  {\rm trace}  \{ \gamma_0 \gamma_5 M^{-1}_d(x,t;0,0)\gamma_0 \gamma_5 M^{-1}_u(0,0;x,t) \}
    $$
\item Euclidean space $\propto A e^{- m t} $
\item Tune bare mass until interacting meson mass is correct, prefactor gives pion, kaon, B meson decay constant
\item etc.. 
\end{itemize}
  
\end{frame}

\begin{frame}[fragile]\small\frametitle{ Hadronic observables}
  \begin{itemize}
    \item Many ways to use Grid to assemble Hadronic observables
    \item Hadrons \link{https://github.com/aportelli/Hadrons}
    \item GPT \link{https://github.com/clehner/gpt}
    \item As a library (e.g. CPS, MILC, Chris Kelly)
    \item Write your own?
  \begin{itemize}
  \item to add {\bf new} functionaly to any of the above, it is necessary to know how to write in Grid
  \item developing simple code is a useful base
  \end{itemize}
  \end{itemize}
\end{frame}

\begin{frame}[fragile]\small\frametitle{ Example code}
  \begin{itemize}
    \item Developed in Grid in about 6h, Mon 30th Aug
    \item Only loosely tested: intended to be illustrative
  \end{itemize}
  \href{https://github.com/paboyle/Grid/blob/develop/examples/Example_Mobius_spectrum.cc}
  {\color{blue}https://github.com/paboyle/Grid/blob/develop/examples/Example\_Mobius\_spectrum.cc}
{\tiny
\begin{verbatim}
  template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
  template<class Field> void GaussianSmear(LatticeGaugeField &U,Field &unsmeared,Field &smeared);
  void PointSource(Coordinate &coor,LatticePropagator &source);
  void Z2WallSource(GridParallelRNG &RNG,int tslice,LatticePropagator &source);
  void GaussianSource(Coordinate &site,LatticeGaugeField &U,LatticePropagator &source);
  void GaussianWallSource(GridParallelRNG &RNG,int tslice,LatticeGaugeField &U,LatticePropagator &source);
  void SequentialSource(int tslice,Coordinate &mom,LatticePropagator &spectator,LatticePropagator &source);
  void MakePhase(Coordinate mom,LatticeComplex &phase);
  template<class Action> void Solve(Action &D,LatticePropagator &source,LatticePropagator &propagator);
  void MesonTrace(std::string file,LatticePropagator &q1,LatticePropagator &q2,LatticeComplex &phase)
\end{verbatim}
}  

\end{frame}

\begin{frame}[fragile]\small\frametitle{ Gist of the programme }
{\tiny
\begin{verbatim}
  LatticePropagator point_source(UGrid);
  LatticePropagator wall_source(UGrid);
  LatticePropagator gaussian_source(UGrid);

  Coordinate Origin({0,0,0,0});
  PointSource   (Origin,point_source);
  Z2WallSource  (RNG4,0,wall_source);
  GaussianSource(Origin,Umu,gaussian_source);
  
  std::vector<LatticePropagator> PointProps(nmass,UGrid);
  std::vector<LatticePropagator> GaussProps(nmass,UGrid);
  std::vector<LatticePropagator> Z2Props   (nmass,UGrid);

  for(int m=0;m<nmass;m++) {
    
    Solve(*FermActs[m],point_source   ,PointProps[m]);
    Solve(*FermActs[m],gaussian_source,GaussProps[m]);
    Solve(*FermActs[m],wall_source    ,Z2Props[m]);
  
  }

  LatticeComplex phase(UGrid);
  Coordinate mom({0,0,0,0});
  MakePhase(mom,phase);
  
  for(int m1=0 ;m1<nmass;m1++) {
  for(int m2=m1;m2<nmass;m2++) {
    std::stringstream ssp,ssg,ssz;

    ssp<<config<< "_m" << m1 << "_m"<< m2 << "_point_meson.xml";
    ssg<<config<< "_m" << m1 << "_m"<< m2 << "_smeared_meson.xml";
    ssz<<config<< "_m" << m1 << "_m"<< m2 << "_wall_meson.xml";

    MesonTrace(ssp.str(),PointProps[m1],PointProps[m2],phase);
    MesonTrace(ssg.str(),GaussProps[m1],GaussProps[m2],phase);
    MesonTrace(ssz.str(),Z2Props[m1],Z2Props[m2],phase);
  }}
\end{verbatim}
}  
\end{frame}

\begin{frame}[fragile]\small\frametitle{Loading a configuration}
  
{\tiny
\begin{verbatim}
  LatticeGaugeField Umu(UGrid);
  std::string config;
  if( argc > 1 && argv[1][0] != '-' )
  {
    std::cout<<GridLogMessage <<"Loading configuration from "<<argv[1]<<std::endl;
    FieldMetaData header;
    NerscIO::readConfiguration(Umu, header, argv[1]);
    config=argv[1];
  }
  else
  {
    std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
    SU<Nc>::ColdConfiguration(Umu);
    //    SU<Nc>::HotConfiguration(RNG4,Umu);
    config="ColdConfig";
  }
\end{verbatim}
}
\begin{itemize}
  \item SciDAC/ILDG and Binary formats too
  \item MPI2IO is used
  \end{itemize}
\end{frame}

\begin{frame}[fragile]\small\frametitle{Fermion action}

{\tiny
\begin{verbatim}
  std::vector<RealD> masses({ 0.03,0.04,0.45} ); // u/d, s, c ??

  int nmass = masses.size();

  std::vector<MobiusFermionR *> FermActs;
  
  std::cout<<GridLogMessage <<"======================"<<std::endl;
  std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl;
  std::cout<<GridLogMessage <<"======================"<<std::endl;

  for(auto mass: masses) {

    RealD M5=1.0;
    RealD b=1.5;// Scale factor b+c=2, b-c=1
    RealD c=0.5;
    
    FermActs.push_back(new MobiusFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c));
   
  }
\end{verbatim}
}
\begin{itemize}
\item Wilson, Clover, Twisted mass, Staggered, other chiral fermions similar
\end{itemize}
\end{frame}

\begin{frame}[fragile]\small\frametitle{Red Black Solver}

{\tiny
\begin{verbatim}
template<class Action>
void Solve(Action &D,LatticePropagator &source,LatticePropagator &propagator)
{
  GridBase *UGrid = D.GaugeGrid();
  GridBase *FGrid = D.FermionGrid();

  LatticeFermion src4  (UGrid); 
  LatticeFermion src5  (FGrid); 
  LatticeFermion result5(FGrid);
  LatticeFermion result4(UGrid);
  
  ConjugateGradient<LatticeFermion> CG(1.0e-8,100000);
  SchurRedBlackDiagMooeeSolve<LatticeFermion> schur(CG);
  ZeroGuesser<LatticeFermion> ZG; // Could be a DeflatedGuesser if have eigenvectors
  for(int s=0;s<Nd;s++){
    for(int c=0;c<Nc;c++){
      PropToFerm<Action>(src4,source,s,c);

      D.ImportPhysicalFermionSource(src4,src5);

      result5=Zero();
      schur(D,src5,result5,ZG);
      std::cout<<GridLogMessage
	       <<"spin "<<s<<" color "<<c
	       <<" norm2(src5d) "   <<norm2(src5)
               <<" norm2(result5d) "<<norm2(result5)<<std::endl;

      D.ExportPhysicalFermionSolution(result5,result4);

      FermToProp<Action>(propagator,result4,s,c);
    }
  }
}
\end{verbatim}
}  
\begin{itemize}
\item Interface works for both 4d and 5d fermion types
\end{itemize}
\end{frame}

\begin{frame}[fragile]\small\frametitle{ Sources  }

{\tiny
\begin{verbatim}
void PointSource(Coordinate &coor,LatticePropagator &source)
{
  source=Zero();
  SpinColourMatrix kronecker; kronecker=1.0;
  pokeSite(kronecker,source,coor);
}
void Z2WallSource(GridParallelRNG &RNG,int tslice,LatticePropagator &source)
{
  GridBase *grid = source.Grid();
  LatticeComplex noise(grid);
  LatticeComplex zz(grid); zz=Zero();
  LatticeInteger t(grid);

  RealD nrm=1.0/sqrt(2);
  bernoulli(RNG, noise); // 0,1 50:50

  noise = (2.*noise - Complex(1,1))*nrm;

  LatticeCoordinate(t,Tdir);
  noise = where(t==Integer(tslice), noise, zz);

  source = 1.0;
  source = source*noise;
}
\end{verbatim}
}  

\end{frame}


\begin{frame}[fragile]\small\frametitle{ Smearing }
  \begin{itemize}
    \item Reuse the smearing we developed in earlier lectures!
  \end{itemize}

{\tiny
\begin{verbatim}
template<class Field>
void GaussianSmear(LatticeGaugeField &U,Field &unsmeared,Field &smeared)
{
  typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
  Laplacian_t Laplacian(U);

  Integer Iterations = 40;
  Real width = 2.0;
  Real coeff = (width*width) / Real(4*Iterations);

  Field tmp(U.Grid());
  smeared=unsmeared;
  //  chi = (1-p^2/2N)^N kronecker
  for(int n = 0; n < Iterations; ++n) {
    Laplacian.M(smeared,tmp);
    smeared = smeared - coeff*tmp;
  }
}
void GaussianSource(Coordinate &site,LatticeGaugeField &U,LatticePropagator &source)
{
  LatticePropagator tmp(source.Grid());
  PointSource(site,source);
  std::cout << " GaussianSource Kronecker "<< norm2(source)<<std::endl;
  tmp = source;
  GaussianSmear(U,tmp,source);
  std::cout << " GaussianSource Smeared "<< norm2(source)<<std::endl;
}
void GaussianWallSource(GridParallelRNG &RNG,int tslice,LatticeGaugeField &U,LatticePropagator &source)
{
  Z2WallSource(RNG,tslice,source);
  auto tmp = source;
  GaussianSmear(U,tmp,source);
}
\end{verbatim}
}  
\end{frame}

\begin{frame}[fragile]\small\frametitle{ Meson three point functions }

  \begin{itemize}
    \item Use sequential source approach; contract an extended propagator via standard meson contraction
    \end{itemize}
{\tiny
\begin{verbatim}
void SequentialSource(int tslice,Coordinate &mom,LatticePropagator &spectator,LatticePropagator &source)
{
  assert(mom.size()==Nd);
  assert(mom[Tdir] == 0);

  GridBase * grid = spectator.Grid();
  Gamma G5(Gamma::Algebra::Gamma5);

  LatticeInteger ts(grid);
  LatticeCoordinate(ts,Tdir);
  source = Zero();
  source = where(ts==Integer(tslice),spectator,source); // Stick in a slice of the spectator, zero everywhere else

  LatticeComplex phase(grid);
  MakePhase(mom,phase);

  source = G5*source *phase;
}
\end{verbatim}
}  

\end{frame}

\begin{frame}[fragile]\small\frametitle{ Contractions }
{\miniscule
\begin{verbatim}
class MesonFile: Serializable {
public:
  GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFile, std::vector<std::vector<Complex> >, data);
};

void MesonTrace(std::string file,LatticePropagator &q1,LatticePropagator &q2,LatticeComplex &phase)
{
  const int nchannel=4;
  Gamma::Algebra Gammas[nchannel][2] = {
    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::Gamma5},
    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::GammaTGamma5},
    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::Gamma5},
    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::GammaTGamma5}
  };

  Gamma G5(Gamma::Algebra::Gamma5);

  LatticeComplex meson_CF(q1.Grid());
  MesonFile MF;

  for(int ch=0;ch<nchannel;ch++){

    Gamma Gsrc(Gammas[ch][0]);
    Gamma Gsnk(Gammas[ch][1]);

    meson_CF = trace(G5*adj(q1)*G5*Gsnk*q2*adj(Gsrc));

    std::vector<TComplex> meson_T;
    sliceSum(meson_CF,meson_T, Tdir);

    int nt=meson_T.size();

    std::vector<Complex> corr(nt);
    for(int t=0;t<nt;t++){
      corr[t] = TensorRemove(meson_T[t]); // Yes this is ugly, not figured a work around
      std::cout << " channel "<<ch<<" t "<<t<<" " <<corr[t]<<std::endl;
    }
    MF.data.push_back(corr);
  }

  {
    XmlWriter WR(file);
    write(WR,"MesonFile",MF);
  }
}
\end{verbatim}
}  
\end{frame}

\begin{frame}[fragile]\small\frametitle{ Hadrons}

  \begin{itemize}
  \item Connects many of these ideas and more in reusable modules
  \item Connect outputs to inputs in dataflow style programming/graphs
  \end{itemize}
\begin{center} \link{https://github.com/aportelli/Hadrons/tree/develop/Hadrons/Modules}\end{center}
\includegraphics[width=0.5\textwidth]{HadronsModules.pdf}
\end{frame}

\begin{frame}[fragile]\small\frametitle{ Final words }
  \begin{itemize}
  \item Aims
  \begin{itemize}
  \item convince you that LQCD software can be elegant, portable and fast
  \item convince you that algorithms can be easy to implement
  \item convince you that code can be elegant, portable and fast
  \item convince you to get your hands dirty!
  \item draw connections between a sample of the core algorithms \& methods of LQCD
  \item keep the exposition simple while still covering the depth
  \end{itemize}
  \item Please provide feedback: what worked and what didn't
  \item I hope you enjoyed the course
  \end{itemize}
\end{frame}


\end{document}