 \documentclass[authoryear,3p,11pt]{elsarticle}

\usepackage{amssymb,amsmath,amsfonts,amstext}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{mathtools}
\usepackage{stmaryrd}
\usepackage{color}
\usepackage{xcolor}
\usepackage{hyperref}
\usepackage{float}
\usepackage{rotating}
\usepackage{bigints}
\newcommand{\X}{\mathbf{X}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\Pb}{\mathbf{P}}
\newcommand{\p}{\boldsymbol \varphi}
\newcommand{\F}{\mathbf{F}}
\newcommand{\N}{\mathbf{N}}
\newcommand{\M}{\text{M}}
\newcommand{\m}{\text{m}}
\usepackage{indentfirst}
\usepackage{rotating}
\usepackage{multirow}
\setlength{\parindent}{2em}

\setlength{\parindent}{4em}
\setlength{\parskip}{0.2cm}

%\renewcommand\includegraphics[2][a]{FIGURE}


\journal{Journal of the Mechanics and Physics of Solids}
\graphicspath {
    {./}
    {./Figures/}
}

\setlength{\parindent}{0pt}
\begin{document} 

\begin{center}README FOR CODE USED IN \\
\textit{Harnessing fluctuations to discover dissipative evolution equations}
\end{center}

\begin{center}{\small Xiaoguai Li, Nicolas Dirr, Peter Embacher, Johannes Zimmer, and Celia Reina}
\end{center}

\begin{center}{\small April 2019}
\end{center}

\vspace{1cm}

\tableofcontents
\newpage

\section{Overview} 

This file explains all codes used in

\emph{Harnessing fluctuations to discover dissipative evolution equations
}, by Xiaoguai Li, Nicolas Dirr, Peter Embacher, Johannes Zimmer, and Celia Reina.

to generate Figures 3-6 as well as the embedded movies. Specifically, the generation of the figures required four steps: simulations of the zero range process, data post-processing to convert particle fluctuations to components of the discretized dissipative operator, fit of these components (independently or accounting for mass conservation as a constraint), and macroscopic simulations. The names of the folders and their corresponding purpose are listed here.

\begin{itemize}
\item {\bf 1\_ZeroRangeProcess}: C++ code to carry out multiple realizations of a zero range process (see Section \ref{Sec:ZRP} for further details).
\item {\bf 2\_Postprocessing}: Matlab code to compute each component of the discretized dissipative operator from particle data for discrete values of $\rho$ and $\nabla \rho$ (see Section \ref{Sec:Postprocessing} for further details).
\item {\bf 3a\_IndependentFit}: Matlab code to obtain a polynomial fit of each independent nonzero component of the discretized operator in the space ($\rho, \nabla \rho$) (see Section \ref{Sec:IndepFit} for further details).
\item {\bf 3b\_FitWithConstraint}: Matlab code to perform a fit to the nonzero components of the discretized operator that satisfies mass conservation --- this is imposed as a constraint --- (see Section \ref{Sec:FitConstraint} for further details).
\item {\bf 4\_MacroSimulations}: Matlab code to perform coarse graining and compute particle-based solution. (Fig.5, Fig.6(b) and supplementary movies) (see Section \ref{Sec:MacroSimulations} for further details).
\end{itemize}

g++11 compiler is needed for the codes in folder 1\_ZeroRangeProcess, and all Matlab codes included in this project were tested with Matlab R2014b.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Zero Range Process (folder 1\_ZeroRangeProcess)} \label{Sec:ZRP}

\subsection{Overview of the code: Main file and functions} 


The folder {\bf 1\_ZeroRangeProcess} is composed of the following files.
\begin{itemize}
\item {\bf Main.cpp}: Main file to be launched where parameters are set and options are selected (see next subsection for details).
\item {\bf SavingParameters.h/.cpp}: Function that saves the parameters set in the Main file.
\item {\bf InitialProfile.h/.cpp}: Function that generates the initial density profile according to the parameters and options selected.
\item {\bf ZRP$\_$KMC$\_$fixedBoundary.h/cpp}: Functions that advances the density profile using a Lattice Kinetic Monte Carlo algorithm, and considering Dirichlet (fixed) boundary conditions.
\item {\bf Postprocessing.h/cpp}: Function that evaluates and saves $\langle \rho_{\epsilon}(t_0)$, $\gamma_1 \rangle$, $\langle \rho_{\epsilon}(t_0+h),\gamma_1 \rangle$, $\langle \rho_{\epsilon}(t_0),\gamma_2 \rangle, \langle \rho_{\epsilon}(t_0+h),\gamma_2 \rangle, \ldots$ for each realization.
\end{itemize}


\subsection{Parameters and options of the code} 


The following variables and options can be adjusted in the Main file. Their equivalent name in the paper and their meaning is provided below.

\begin{itemize}

\item {\bf Configuration and process parameters}:\\
\begin{tabular}{ p{2.2cm} p{12.5cm} }
  process & this option can only be set to ``ZRP''. \\
  profile & this option can be set as ``flat'', ``linear'' or ``triangle'' to generate a flat profile, a linear profile or a triangular profile. Note that the profile is randomly perturbed in InitialProfile.cpp. Flat profile uses parameters ``Nbin'' and ``Npart''. ``Linear'' and ``triangle'' profile uses parameters ``Nbin'', ``slope'' and ``center''.\\
    Nbin & number of bins  \\
  Npart & number of particles  \\
  slope & this is only useful when profile is set to be ``linear'' or ``triangle''. It sets slope of linear profile, or the positive slope for the triangular profile (the magnitudes of positive and negative slopes are identical in the triangular profile).\\
center & this is only useful when profile is set to be ``linear" or ``triangle". It is the value at the center of the profile.
\end{tabular}

\item {\bf Macroscopic time parameters}: \\
\begin{tabular}{p{2.2cm} p{12.5cm}  }
  t$\_$equilibration &  $t_{prep}-t_{ini}$, time of equilibration. \\
  t$\_$randomize &   $t_0 - t_{prep}$, time to prepare the system for a new realization.\\
  h & $h$, simulation time for actual measurements.\\
\end{tabular}

\item {\bf Sampling parameters}:\\
\begin{tabular}{p{2.2cm} p{12.5cm}  }
 R1 &  $R_1$ \\
 R2 &   $R_2$\\
\end{tabular}

\item {\bf Basis function parameters}:\\
\begin{tabular}{p{2.2cm} p{12.5cm}  }
Ngamma1 & Number of gamma functions in the first basis.\\
Ngamma2 & Number of gamma functions in the second basis.\\
&  Note: The code is set to postprocess the same ZRP data with two sets of basis functions, although only results from one of them is shown in the paper.
\end{tabular}

\item {\bf Output options}:\\
\begin{tabular}{p{2.2cm} p{12.5cm}  }
 file$\_$path1 &  Directory where the user wants the save the output of the first group of basis functions. \\
 file$\_$path2 &  Directory where the user wants the save the output of the second group of basis functions.
\end{tabular}
\end{itemize}

In addition, the type of Zero Range Process can be set inside the function ZRP$\_$KMC$\_$fixedBoundary.cpp, by adjusting the parameter ``power'' to 1.0 or 2.0, for instance, for g(k)=k and g(k)=k$^2$, respectively. All simulations in the article used the value 2.0.


\subsection{Running the code} 


The following commands may be typed on the terminal to compile the code

{\tt >> module load gcc/6.3.0 } \\
{\tt >> g++ -std=gnu++0x -O3 Main.cpp InitialProfile.cpp ZRP$\_$KMC$\_$new.cpp Postprocessing.cpp SavingParameters.cpp -o executable$\_$name}

where the first command is needed to load version C++11 in order to use the random variable generator used in the code.


\subsection{Output of the code and postprocessing} \label{Sec:OutputZRP}


After running the executable, the code generates in each filename path, one summary file and R1 data files:

\begin{itemize}
\item {\bf Summary.m}: contains all the parameters prescribed in the Main file.
\item {\bf Data$\_0$.m}: each row contains $\langle \rho_{\epsilon}(t_0)$, $\gamma_1 \rangle$, $\langle \rho_{\epsilon}(t_0+h),\gamma_1 \rangle$, $\langle \rho_{\epsilon}(t_0),\gamma_2 \rangle, \langle \rho_{\epsilon}(t_0+h),\gamma_2 \rangle, \ldots$ for each realization $r\in[1:$R$2]$, i.e.~there are R2 rows. The number in this data file name ranges from $0$ to R$1$.
\end{itemize}

This data compiles the information of all realizations of a single profile, and can be further postprocessed using the code in folder {\bf 2$\_$Postprocessing} to deliver the expected values of $\rho$ and $\nabla \rho$, and each nonzero component of the dissipative operator discretized using the $\gamma$ functions. We remark that each profile must be postprocessed independently.

\subsection{Parameters used to generate the figures in the paper}

The default parameters for Figures 4-6 are:
\begin{center}
\begin{tabular}{  | l |c|c|c|c|l|c|c|c|}
\hline
Nbin &t\_equilibration &t\_randomize &h &R1 &R2 &Ngamma &x\_basis \\
\hline
5000	&4e-06	&4e-09	&4e-11	&400	&2000	&40	&[0.0:0.025:0.975] \\
\hline
\end{tabular}
\end{center}

Profile: ``flat'' was used with Npart = [20000:500:50000]; For the data with conservation constraint, ``linear'' was used with slope = $\pm$[5, 11] and center = 7,  and ``triangle''  was used with slope = $\pm$[15, 19] and center = 7; For the data without conservation constraint, ``linear'' was used with slope = [-13:1:13] and center = 7, and ``triangle'' was used with slope = [14:1:19] and [-19:1:-14], and center = 7.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Postprocessing (folder 2\_Postprocessing)} \label{Sec:Postprocessing}

\subsection{Code and test data} \label{Sec:PostCode}
Folder {\bf 2\_Postprocessing} includes the following files.
\begin{itemize}
\item {\bf Postprocessing.m}: the Matlab script used to do postprocessing of the data from the zero range process, and estimate the elements of the operator that are expected to be non-zero.
%compute the expected value of each nonzero component of the discretized operator. 
\item {\bf Summary.m}: output file from the previous code, cf.~Section \ref{Sec:OutputZRP}. $\bf Postprocessing.m$ reads this file, and writes new data at the end of it after postprocessing. Specifically, it provides the expected value for the following quantities for each pair of neighboring $\gamma$ functions: 
 \begin{itemize}
 \item[] $\rho_b$, $\nabla \rho_b$, $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_b \rangle$, using variable names {\bf rho\_b\_b},  {\bf drho\_b\_b} and {\bf diagonal\_b\_b}, respectively.
    \item[] $\rho_b$, $\nabla \rho_b$, $\langle \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)} \gamma_b,\gamma_{b+1} \rangle$, using variable names {\bf rho\_b\_bPlus1}, {\bf drho\_b\_bPlus1} and {\bf off\_b\_bPlus1}
  \item[] $\rho_b$, $\nabla \rho_b$, $\langle \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)} \gamma_b,\gamma_{b-1} \rangle$, using variable names {\bf rho\_b\_bMinus1}, {\bf drho\_b\_bMinus1} and {\bf off\_b\_bMinus1}
 \end{itemize}
 We remark that for a given pair of neighboring $\gamma$ functions, these are denoted as $\gamma_b$ and $\gamma_{b+1}$ (from left to right), for the purposes of computing $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_b \rangle$ and \\ $\langle \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)} \gamma_b,\gamma_{b+1} \rangle$, and the same functions are labeled as $\gamma_{b-1}$ and $\gamma_b$ (from left to right), for the purpose of computing $\langle \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)} \gamma_b,\gamma_{b-1} \rangle$. For this reason, the vectors {\bf rho\_b\_b}, {\bf rho\_b\_bPlus1} and {\bf rho\_b\_bMinus1} do not coincide.
\item {\bf Data\_0.m, ..., Data\_4.m}: Output files from the previous code, cf.~Section \ref{Sec:OutputZRP}. 
\end{itemize}

The files {\bf Summary.m}, {\bf Data\_0.m, ..., Data\_4.m} provided in the folder correspond to an example with $R_1=5$ and $R_2=200$, which may be used to test the code {\bf Postprocessing.m}. The real dataset is too large to be attached here. 

\subsection{Output of the code and further steps} 
The output is added to $\bf Summary.m$, as previously explained.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Independent fit of the components (folder 3a\_IndependentFit)} \label{Sec:IndepFit}

\subsection{Real dataset from the zero range process after postprocessing} \label{Sec:IndepFitData}

The previous code delivers a {\bf Summary.m} file for a specific density profile, that contains the three nonzero components of the discretized operators as a function of discrete values $\rho$ and $\nabla \rho$, probed with such profile. The results from the various profiles must then be compiled into single vectors {\bf rho\_b\_b},  {\bf drho\_b\_b}, {\bf diagonal\_b\_b}, {\bf rho\_b\_bPlus1}, {\bf drho\_b\_bPlus1}, {\bf off\_b\_bPlus1}  {\bf rho\_b\_bMinus1}, {\bf drho\_b\_bMinus1}, {\bf off\_b\_bMinus1} and saved in a  $\bf Data.mat$  Matlab file. The file $\bf Data.mat$ included in folder {\bf 3a\_IndependentFit} corresponds to the actual data used in the paper.

\subsection{Code overview} 
There are three scripts used to fit three nonzero components in discretized operator independently. Their names and purposes are listed below.
\begin{itemize}
\item {\bf fit\_dia\_b\_b.m}: polynomial fit of $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_b \rangle$ as a function of $\rho_b$ and $\nabla \rho_b$.
\item {\bf fit\_off\_b\_bPlus1.m}: polynomial fit of $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_{b+1} \rangle$ as a function of  $\rho_b$ and $\nabla \rho_b$.
\item {\bf fit\_off\_b\_bMinus1.m}: polynomial fit of $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_{b-1} \rangle$ as a function of $\rho_b$ and $\nabla \rho_b$.
\end{itemize}

\subsection{Options of the code} 
There are three fit options in these three scripts, which can be set by parameter 'fitID', and determines the order of the polynomial in the variables $\rho_b$ and $\nabla \rho_b$.

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
fitID  & order in $\rho_b$ &  order in $\nabla \rho_b$ \\
\hline
 22                & 2          & 2\\
\hline
21	& 2          & 1\\
\hline
20	& 2          & 0\\
\hline
\end{tabular}
\end{center}

In Figure 4, all nonzero components are fit with fitID=22. 

\subsection{Output of the code} 
These three scripts give the following outputs. 

\begin{center}
\begin{tabular}{ |l|l|p{76mm}| }
\hline
\multirow{3}{*}{fit parameters} &  {\bf dia\_fit\_22.mat} & {\centering Fit parameters of $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_{b} \rangle$ with fitID=22. }\\ 
\cline{2-3}
& {\bf off\_b\_bMinus1\_fit\_22.mat} &  Fit parameters of  $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_{b-1} \rangle$ with fitID=22. \\ 
\cline{2-3}
& {\bf off\_b\_bPlus1\_fit\_22.mat} &  Fit parameters of  $\langle  \mathcal{K}_{\left(\rho_b+ \nabla \rho|_b (x-x_b) +\ldots\right)}\gamma_b,\gamma_{b+1} \rangle$ with fitID=22. \\ 
\hline
\multirow{3}{*}{{\bf figures}} &  {\bf dia\_b\_b.fig} & Plots of the raw data and fitted surface. \\ 
& {\bf off\_b\_bMinus1.fig} &  \\ 
& {\bf off\_b\_bPlus1.fig} & \\ 
\hline
\multirow{2}{*}{{\bf errorPlot}} &  {\bf error\_dia\_22.fig} & \multirow{3}{75mm}{Plots of the relative error between the fitted results and analytic values, known for the zero range process considered. In practice, analytic values of the nonzero components are not necessarily known.} \\
& {\bf error\_off\_b\_bMinus1\_22.fig}  & \\ 
& {\bf error\_off\_b\_bPlus1\_22.fig} &  \\ [10mm] 
\hline
\end{tabular}
\end{center}

Files {\bf dia\_fit\_22.mat}, {\bf off\_b\_bMinus1\_fit\_22.mat} and {\bf off\_b\_bPlus1\_fit\_22.mat} will be used to carry out the macroscopic simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Fit of components with mass conservation constraint (folder 3b\_FitWithConstraint)} \label{Sec:FitConstraint}

\subsection{Real dataset from zero range process after postprocessing} 

The file $\bf Data.mat$ included in folder {\bf 3b\_FitWithConstraint} contains the actual data used in Fig.~6 of the paper using the following scripts, and corresponds to a subset of the data included in  {\bf3a\_IndependentFit/Data.mat}. The variable names are identical to the ones of the previous Section, and these are explained in Subsections \ref{Sec:PostCode} and \ref{Sec:IndepFitData}. 

\subsection{Code and parameters} 
There is a single script in this folder, $\bf optimization\_threeComponents.m$, which is the Matlab script used to fit the three nonzero components together so as to guarantee mass conservation. Specifically, this conservation constraint implies that sum of these three components is zero, which is directly imposed by reducing the number of unknown polynomial coefficients (similarly to before, a quadratic polynomial in both $\rho_b$ and $\nabla \rho_b$ is considered). An objective error function is defined as the $L_2$ norm of the difference between the raw data and the quadratic surface, and this is minimized using a quasi-Newton method to identify the optimal polynomial coefficients.

\subsection{Output of the code} 
These three scripts give the following outputs. 

\begin{center}
\begin{tabular}{ |l|p{105mm}| }
\hline
{\bf fitParameter.mat} & Fit parameters for the quadratic polynomial target functions.\\
\hline
{\bf error\_dia\_22.fig} & \multirow{3}{105mm}{Plots of the relative error between the fitted results and analytic values, known for the zero range process considered. In practice, analytic values of the nonzero components are not necessarily known.} \\
{\bf error\_off\_b\_bMinus1\_22.fig}  & \\ 
{\bf error\_off\_b\_bPlus1\_22.fig} &  \\[6mm]
\hline
\end{tabular}
\end{center}

The file {\bf fitParameter.mat} will serve as input data for the macroscopic simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Macroscopic simulations (folder 4\_MacroSimulations)}  \label{Sec:MacroSimulations}
\subsection{Code overview} 
In folder {\bf 4\_MacroSimulations} there are two Matlab scripts used to perform coarse graining. Their names and purposes are listed below.

\begin{itemize}
\item {\bf MacroEvolution.m}: main script that delivers the macroscopic evolution using the particle-based operator and the analytical operator.
\item {\bf prepare\_analytic\_m.m}: function used in main script, which prepares the analytic mobility $m$ used in the computation of the analytic solution.
\item {\bf grad.m}: function used in main script, which computes gradient for an array, with consideration of periodic boundary condition.
\end{itemize}

\subsection{Options of the code} 
There are two main options in the code that enable the user to choose between the constrained or unconstrained fit of the discrete operator components, and to specify the type of boundary conditions to be used in the macroscopic simulations. In all cases, a macroscopic domain $[0,1]$ is considered in Figs.~5 and 6 of the article.

\begin{center}
\begin{tabular}{ |l|l|p{90mm}| }
\hline
\multirow{3}{*}{constraint} &  0 & Reads fit result from independent fit.\\ 
\cline{2-3}
& 1 &  Reads fit result from fit with conservation constraint.\\  
\hline
\multirow{3}{*}{boundaryCondition} &  'periodic' & Periodic boundary conditions.\\ 
\cline{2-3}
& 'Dirichlet' & Dirichlet boundary conditions.\\  
\hline
\end{tabular}
\end{center}

\subsection{Parameters of the code} 

\begin{itemize}
\item {\bf Discretization parameters}:

\begin{tabular}{ p{2.2cm} p{12.5cm} }
Ngamma &  number of basis functions. This value needs to be in agreement with the value of Ngamma used to postprocess the particle data.\\
dx & mesh size (dx=1/Ngamma)\\
dt & discretization in time\\
T & total simulation time\\
timesteps &  total simulation steps (timesteps=T/dt)\\
%x & position of discretized points\\
%N & number of discretized points
\end{tabular}

\item {\bf Initial profile}: 

\begin{tabular}{p{2.2cm} p{12.5cm}  }
density\_0 &  initial profile (t=0) as a function of x 
\end{tabular}

\item {\bf Video setting}:

\begin{tabular}{p{2.2cm} p{12.5cm}  }
dump &  number of frames\\
v &   video name and format\\
v.FrameRate & frame speed
\end{tabular}

\item {\bf Output snapshots}:

\begin{tabular}{p{2.2cm} p{12.5cm}  }
snapshots & specify time to output snapshots; real simulation time is set to be snapshots/1000
\end{tabular}

\end{itemize}

\subsection{Output of the code} 
Output of the code includes two parts:
\begin{itemize}
\item {\bf movies}: this code automatically makes a directory named 'movies', and saves in it a movie of the analytical evolution and the particle-based evolution.
\item {\bf snapshots}: this code automatically makes a directory named 'snapshots', and saves in it the data of the analytic profile and the particle-based profile at the inquired times. An independent data file is created for each inquired time, and their names coincide with the values entered in the vector snapshots (see parameters of the code).
\end{itemize}

\subsection{Options and parameters used to generate the figures in the paper}

$\bf Options:$
\begin{center}
\begin{tabular}{ |c|c|c|c|c| }
\hline
Figures; Movies & constraint  & boundaryCondition  & snapshots & density\_0 \\
\hline
Fig.~5(a); Movie1 &  0 &  'periodic' & [1 3 6 10] &-3*cos(2*pi*x/xfinal)+7\\
\hline
Fig.~5(b); Movie2 & 0 &'Dirichlet' & [1 3 8] &-5*cos(3*pi*x/xfinal)+7\\
\hline
Fig.~6(b); Movie3 & 1 &'periodic' & [1 3 6 10]&-3*cos(2*pi*x/xfinal)+7\\
\hline
\end{tabular}
\end{center}


$\bf Common Parameters:$
\begin{center}
\begin{tabular}{ |c|c|c|c|c|c| }
\hline
Ngamma & dx & dt & T & dump & v.FrameRate\\
\hline
40 & 0.025 & 6.2500e-07 & 0.01 &100 & 5 \\
\hline
\end{tabular}
\end{center}

%\vskip 0.8in
%\begin{center}
%\textit{Enjoy the code!}
%\end{center}

\end{document}

\endinput 
