Ian Jauslin

Ian Jauslin

Home Publications Talks Software Teaching Contact
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile54
-rw-r--r--docs/nstrophy_doc/Makefile70
-rw-r--r--docs/nstrophy_doc/bibliography/conf.BBlog5
-rw-r--r--docs/nstrophy_doc/libs/ian.cls659
-rw-r--r--docs/nstrophy_doc/libs/iantheo.sty162
-rw-r--r--docs/nstrophy_doc/libs/largearray.sty19
-rw-r--r--docs/nstrophy_doc/libs/point.sty106
-rw-r--r--docs/nstrophy_doc/nstrophy_doc.tex146
-rw-r--r--src/main.c282
-rw-r--r--src/navier-stokes.c299
-rw-r--r--src/navier-stokes.h48
11 files changed, 1850 insertions, 0 deletions
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..36a266c
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,54 @@
+PROJECTNAME=nstrophy
+VERSION=0.0
+
+#DB= -ggdb
+OPT= -O3
+
+WARNINGS=-Wall -Wextra -Wno-strict-overflow -std=c99 -Wpedantic
+
+PREFIX=/usr
+BINDIR=$(PREFIX)/bin
+MANDIR=$(PREFIX)/share/man/man1
+
+CC=/usr/bin/gcc
+LD=$(CC)
+
+#INCLUDES =
+
+#LIBDIRS =
+LIBS = -lm -lfftw3
+
+
+override LDFLAGS +=$(LIBDIRS)$(LIBS)
+override CFLAGS +=$(INCLUDES)$(DB)
+
+override CFLAGS +=$(OPT) $(WARNINGS)
+
+BUILDDIR=./build
+SRCDIR=./src
+
+OBJS=$(patsubst %.c,%.o,$(wildcard $(SRCDIR)/*.c))
+
+all: dist
+
+$(PROJECTNAME): $(OBJS)
+ $(LD) -o $@ $^ $(LDFLAGS)
+
+%.o: %.c
+ $(CC) -c $(CFLAGS) $< -o $@
+
+dist: $(PROJECTNAME)
+ /bin/rm -rf $(BUILDDIR)
+ /bin/mkdir $(BUILDDIR)
+ /bin/cp $(PROJECTNAME) $(BUILDDIR)
+ @rm -f $(SRCDIR)/*.o
+ @rm -f $(PROJECTNAME)
+
+install: dist
+ install -Dm755 $(BUILDDIR)/$(PROJECTNAME) $(BINDIR)/$(PROJECTNAME)
+ install -Dm644 man/$(PROJECTNAME).1.gz $(MANDIR)/$(PROJECTNAME).1.gz
+
+clean:
+ @rm -f $(SRCDIR)/*.o
+ @rm -rf $(BUILDDIR)
+ @rm -f $(PROJECTNAME)
diff --git a/docs/nstrophy_doc/Makefile b/docs/nstrophy_doc/Makefile
new file mode 100644
index 0000000..c63e53f
--- /dev/null
+++ b/docs/nstrophy_doc/Makefile
@@ -0,0 +1,70 @@
+PROJECTNAME=$(basename $(wildcard *.tex))
+LIBS=$(notdir $(wildcard libs/*))
+FIGS=$(notdir $(wildcard figs/*.fig))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+SYNCTEXS=$(addsuffix .synctex.gz, $(PROJECTNAME))
+
+all: $(PROJECTNAME)
+
+$(PROJECTNAME): $(LIBS) $(FIGS)
+ pdflatex -file-line-error $@.tex
+ pdflatex -file-line-error $@.tex
+ pdflatex -synctex=1 $@.tex
+
+$(PROJECTNAME).aux: $(LIBS) $(FIGS)
+ pdflatex -file-line-error -draftmode $(PROJECTNAME).tex
+
+
+$(SYNCTEXS): $(LIBS) $(FIGS)
+ pdflatex -synctex=1 $(patsubst %.synctex.gz, %.tex, $@)
+
+install: $(PROJECTNAME)
+ cp $^.pdf $(INSTALLDIR)/
+
+install-synctex: $(SYNCTEXS)
+ cp $^ $(INSTALLDIR)/
+
+install-bibliography: bibliography/bibliography.tex
+ cp bibliography/bibliography.tex $(INSTALLDIR)/bibliography/
+
+
+libs: $(LIBS)
+
+$(LIBS):
+ ln -fs libs/$@ ./
+
+bibliography/bibliography.tex: $(PROJECTNAME).aux
+ BBlog -c bibliography/conf.BBlog -d $(BIBLIOGRAPHY) -b bibliography/bibliography.tex
+
+figs: $(FIGS)
+
+$(FIGS):
+ make -C figs/$@
+ ln -fs figs/$@/*.pdf ./
+
+
+clean-aux: clean-figs-aux
+ rm -f $(addsuffix .aux, $(PROJECTNAME))
+ rm -f $(addsuffix .log, $(PROJECTNAME))
+ rm -f $(addsuffix .out, $(PROJECTNAME))
+ rm -f $(addsuffix .toc, $(PROJECTNAME))
+
+clean-libs:
+ rm -f $(LIBS)
+
+clean-figs:
+ $(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean; )
+ rm -f $(notdir $(wildcard figs/*.fig/*.pdf))
+
+clean-figs-aux:
+ $(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean-aux; )
+
+
+clean-tex:
+ rm -f $(PDFS) $(SYNCTEXS)
+
+clean-bibliography:
+ rm -f bibliography/bibliography.tex
+
+clean: clean-aux clean-tex clean-libs clean-figs
diff --git a/docs/nstrophy_doc/bibliography/conf.BBlog b/docs/nstrophy_doc/bibliography/conf.BBlog
new file mode 100644
index 0000000..d1c22bc
--- /dev/null
+++ b/docs/nstrophy_doc/bibliography/conf.BBlog
@@ -0,0 +1,5 @@
+format: \bibitem[%token%]{%citeref%}%author:auth% - {\it %title%}, %journal%, %year%|doi?,\par\penalty10000%n%doi:{\tt\color{blue}\href{http://dx.doi.org/%doi%}{%doi%}}||arxiv?, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/%arxiv%}{%arxiv%}}|.\par\medskip%n% %n%
+out_file: bibliography.tex
+filter:auth: s/([A-Z])[^, ]* /\1. /g; s/ ([^ ,]*),/_\1,_/g; s/ ([^ ,]*)$/_\1/g; s/ //g; s/_/ /g;
+filter:journal: s/([a-zA-Z]) ([0-9]+)/\1~\\-\2/g;
+aux_cmd: \\citation{
diff --git a/docs/nstrophy_doc/libs/ian.cls b/docs/nstrophy_doc/libs/ian.cls
new file mode 100644
index 0000000..538a108
--- /dev/null
+++ b/docs/nstrophy_doc/libs/ian.cls
@@ -0,0 +1,659 @@
+%%
+%% Ian's class file
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% class name
+\ProvidesClass{ian}[2017/09/29]
+
+%% boolean to signal that this class is being used
+\newif\ifianclass
+
+%% options
+% no section numbering in equations
+\DeclareOption{section_in_eq}{\sectionsineqtrue}
+\DeclareOption{section_in_fig}{\sectionsinfigtrue}
+\DeclareOption{section_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{section_in_all}{\sectionsineqtrue\sectionsinfigtrue\PassOptionsToPackage{section_in_theo}{iantheo}}
+\DeclareOption{subsection_in_eq}{\subsectionsineqtrue}
+\DeclareOption{subsection_in_fig}{\subsectionsinfigtrue}
+\DeclareOption{subsection_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{subsection_in_all}{\subsectionsineqtrue\subsectionsinfigtrue\PassOptionsToPackage{subsection_in_theo}{iantheo}}
+\DeclareOption{no_section_in_eq}{\sectionsineqfalse}
+\DeclareOption{no_section_in_fig}{\sectionsinfigfalse}
+\DeclareOption{no_section_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{no_section_in_all}{\sectionsineqfalse\sectionsinfigfalse\PassOptionsToPackage{no_section_in_theo}{iantheo}}
+\DeclareOption{no_subsection_in_eq}{\subsectionsineqfalse}
+\DeclareOption{no_subsection_in_fig}{\subsectionsinfigfalse}
+\DeclareOption{no_subsection_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{no_subsection_in_all}{\subsectionsineqfalse\subsectionsinfigfalse\PassOptionsToPackage{no_subsection_in_theo}{iantheo}}
+% reset point
+\DeclareOption{point_reset_at_section}{\PassOptionsToPackage{reset_at_section}{point}}
+\DeclareOption{point_no_reset_at_section}{\PassOptionsToPackage{no_reset_at_section}{point}}
+\DeclareOption{point_reset_at_theo}{\PassOptionsToPackage{reset_at_theo}{point}}
+\DeclareOption{point_no_reset_at_theo}{\PassOptionsToPackage{no_reset_at_theo}{point}}
+
+\def\ian@defaultoptions{
+ \ExecuteOptions{section_in_all, no_subsection_in_all}
+ \ProcessOptions
+
+ %% required packages
+ \RequirePackage{iantheo}
+ \RequirePackage{point}
+ \RequirePackage{color}
+ \RequirePackage{marginnote}
+ \RequirePackage{amssymb}
+ \PassOptionsToPackage{hidelinks}{hyperref}
+ \RequirePackage{hyperref}
+}
+
+%% paper dimensions
+\setlength\paperheight{297mm}
+\setlength\paperwidth{210mm}
+
+%% fonts
+\input{size11.clo}
+\DeclareOldFontCommand{\rm}{\normalfont\rmfamily}{\mathrm}
+\DeclareOldFontCommand{\sf}{\normalfont\sffamily}{\mathsf}
+\DeclareOldFontCommand{\tt}{\normalfont\ttfamily}{\mathtt}
+\DeclareOldFontCommand{\bf}{\normalfont\bfseries}{\mathbf}
+\DeclareOldFontCommand{\it}{\normalfont\itshape}{\mathit}
+
+%% text dimensions
+\hoffset=-50pt
+\voffset=-72pt
+\textwidth=460pt
+\textheight=704pt
+
+
+%% remove default indentation
+\parindent=0pt
+%% indent command
+\def\indent{\hskip20pt}
+
+%% something is wrong with \thepage, redefine it
+\gdef\thepage{\the\c@page}
+
+%% array lines (to use the array environment)
+\setlength\arraycolsep{5\p@}
+\setlength\arrayrulewidth{.4\p@}
+
+
+%% correct vertical alignment at the end of a document
+\AtEndDocument{
+ \vfill
+ \eject
+}
+
+
+%% hyperlinks
+% hyperlinkcounter
+\newcounter{lncount}
+% hyperref anchor
+\def\hrefanchor{%
+\stepcounter{lncount}%
+\hypertarget{ln.\thelncount}{}%
+}
+
+%% define a command and write it to aux file
+\def\outdef#1#2{%
+ % define command%
+ \expandafter\xdef\csname #1\endcsname{#2}%
+ % hyperlink number%
+ \expandafter\xdef\csname #1@hl\endcsname{\thelncount}%
+ % write command to aux%
+ \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1\endcsname{\csname #1\endcsname}}%
+ \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1@hl\endcsname{\thelncount}}%
+}
+
+%% can call commands even when they are not defined
+\def\safe#1{%
+ \ifdefined#1%
+ #1%
+ \else%
+ {\color{red}\bf?}%
+ \fi%
+}
+
+%% define a label for the latest tag
+%% label defines a command containing the string stored in \tag
+\def\deflabel{
+ \def\label##1{\expandafter\outdef{label@##1}{\safe\tag}}
+
+ \def\ref##1{%
+ % check whether the label is defined (hyperlink runs into errors if this check is omitted)
+ \ifcsname label@##1@hl\endcsname%
+ \hyperlink{ln.\csname label@##1@hl\endcsname}{{\color{blue}\safe\csname label@##1\endcsname}}%
+ \else%
+ \ifcsname label@##1\endcsname%
+ {\color{blue}\csname ##1\endcsname}%
+ \else%
+ {\bf ??}%
+ \fi%
+ \fi%
+ }
+}
+
+
+%% make a custom link at any given location in the document
+\def\makelink#1#2{%
+ \hrefanchor%
+ \outdef{label@#1}{#2}%
+}
+
+
+%% section command
+% counter
+\newcounter{sectioncount}
+% space before section
+\newlength\secskip
+\setlength\secskip{40pt}
+% a prefix to put before the section number, e.g. A for appendices
+\def\sectionprefix{}
+% define some lengths
+\newlength\secnumwidth
+\newlength\sectitlewidth
+\def\section#1{
+ % reset counters
+ \stepcounter{sectioncount}
+ \setcounter{subsectioncount}{0}
+ \ifsectionsineq
+ \setcounter{seqcount}0
+ \fi
+ \ifsectionsinfig
+ \setcounter{figcount}0
+ \fi
+
+ % space before section (if not first)
+ \ifnum\thesectioncount>1
+ \vskip\secskip
+ \penalty-1000
+ \fi
+
+ % hyperref anchor
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\sectionprefix\thesectioncount}
+
+ % get widths
+ \def\@secnum{{\bf\Large\sectionprefix\thesectioncount.\hskip10pt}}
+ \settowidth\secnumwidth{\@secnum}
+ \setlength\sectitlewidth\textwidth
+ \addtolength\sectitlewidth{-\secnumwidth}
+
+ % print name
+ \parbox{\textwidth}{
+ \@secnum
+ \parbox[t]{\sectitlewidth}{\Large\bf #1}}
+
+ % write to table of contents
+ \iftoc
+ % save lncount in aux variable which is written to toc
+ \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@sec.\thesectioncount\endcsname{\thelncount}}
+ \write\tocoutput{\noexpand\tocsection{#1}{\thepage}}
+ \fi
+
+ %space
+ \par\penalty10000
+ \bigskip\penalty10000
+}
+
+%% subsection
+% counter
+\newcounter{subsectioncount}
+% space before subsection
+\newlength\subsecskip
+\setlength\subsecskip{30pt}
+\def\subsection#1{
+ % counters
+ \stepcounter{subsectioncount}
+ \setcounter{subsubsectioncount}{0}
+ \ifsubsectionsineq
+ \setcounter{seqcount}0
+ \fi
+ \ifsubsectionsinfig
+ \setcounter{figcount}0
+ \fi
+
+ % space before subsection (if not first)
+ \ifnum\thesubsectioncount>1
+ \vskip\subsecskip
+ \penalty-500
+ \fi
+
+ % hyperref anchor
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount}
+
+ % get widths
+ \def\@secnum{{\bf\large\hskip.5cm\sectionprefix\thesectioncount.\thesubsectioncount.\hskip5pt}}
+ \settowidth\secnumwidth{\@secnum}
+ \setlength\sectitlewidth\textwidth
+ \addtolength\sectitlewidth{-\secnumwidth}
+ % print name
+ \parbox{\textwidth}{
+ \@secnum
+ \parbox[t]{\sectitlewidth}{\large\bf #1}}
+
+ % write to table of contents
+ \iftoc
+ % save lncount in aux variable which is written to toc
+ \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsec.\thesectioncount.\thesubsectioncount\endcsname{\thelncount}}
+ \write\tocoutput{\noexpand\tocsubsection{#1}{\thepage}}
+ \fi
+
+ % space
+ \par\penalty10000
+ \medskip\penalty10000
+}
+
+%% subsubsection
+% counter
+\newcounter{subsubsectioncount}
+% space before subsubsection
+\newlength\subsubsecskip
+\setlength\subsubsecskip{20pt}
+\def\subsubsection#1{
+ % counters
+ \stepcounter{subsubsectioncount}
+
+ % space before subsubsection (if not first)
+ \ifnum\thesubsubsectioncount>1
+ \vskip\subsubsecskip
+ \penalty-500
+ \fi
+
+ % hyperref anchor
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount}
+
+ % get widths
+ \def\@secnum{{\bf\hskip1.cm\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\hskip5pt}}
+ \settowidth\secnumwidth{\@secnum}
+ \setlength\sectitlewidth\textwidth
+ \addtolength\sectitlewidth{-\secnumwidth}
+ % print name
+ \parbox{\textwidth}{
+ \@secnum
+ \parbox[t]{\sectitlewidth}{\large\bf #1}}
+
+ % write to table of contents
+ \iftoc
+ % save lncount in aux variable which is written to toc
+ \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsubsec.\thesectioncount.\thesubsectioncount.\thesubsubsectioncount\endcsname{\thelncount}}
+ \write\tocoutput{\noexpand\tocsubsubsection{#1}{\thepage}}
+ \fi
+
+ % space
+ \par\penalty10000
+ \medskip\penalty10000
+}
+
+%% itemize
+\newlength\itemizeskip
+% left margin for items
+\setlength\itemizeskip{20pt}
+\newlength\itemizeseparator
+% space between the item symbol and the text
+\setlength\itemizeseparator{5pt}
+% penalty preceding an itemize
+\newcount\itemizepenalty
+\itemizepenalty=0
+% counter counting the itemize level
+\newcounter{itemizecount}
+
+% item symbol
+\def\itemizept#1{
+ \ifnum#1=1
+ \textbullet
+ \else
+ $\scriptstyle\blacktriangleright$
+ \fi
+}
+
+
+\newlength\current@itemizeskip
+\setlength\current@itemizeskip{0pt}
+\def\itemize{%
+ \par\expandafter\penalty\the\itemizepenalty\medskip\expandafter\penalty\the\itemizepenalty%
+ \addtocounter{itemizecount}{1}%
+ \addtolength\current@itemizeskip{\itemizeskip}%
+ \leftskip\current@itemizeskip%
+}
+\def\enditemize{%
+ \addtocounter{itemizecount}{-1}%
+ \addtolength\current@itemizeskip{-\itemizeskip}%
+ \par\expandafter\penalty\the\itemizepenalty\leftskip\current@itemizeskip%
+ \medskip\expandafter\penalty\the\itemizepenalty%
+}
+\newlength\itempt@total
+\def\item{
+ \settowidth\itempt@total{\itemizept\theitemizecount}
+ \addtolength\itempt@total{\itemizeseparator}
+ \par
+ \medskip
+ \hskip-\itempt@total\itemizept\theitemizecount\hskip\itemizeseparator
+}
+
+
+%% prevent page breaks after itemize
+\newcount\previtemizepenalty
+\def\nopagebreakafteritemize{
+ \previtemizepenalty=\itemizepenalty
+ \itemizepenalty=10000
+}
+%% back to previous value
+\def\restorepagebreakafteritemize{
+ \itemizepenalty=\previtemizepenalty
+}
+
+%% enumerate
+\newcounter{enumerate@count}
+\def\enumerate{
+ \setcounter{enumerate@count}0
+ \let\olditem\item
+ \let\olditemizept\itemizept
+ \def\item{
+ % counter
+ \stepcounter{enumerate@count}
+ % set header
+ \def\itemizept{\theenumerate@count.}
+ % hyperref anchor
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\theenumerate@count}
+ \olditem
+ }
+ \itemize
+}
+\def\endenumerate{
+ \enditemize
+ \let\item\olditem
+ \let\itemizept\olditemizept
+}
+
+
+%% equation numbering
+% counter
+\newcounter{seqcount}
+% booleans (write section or subsection in equation number)
+\newif\ifsectionsineq
+\newif\ifsubsectionsineq
+\def\seqcount{
+ \stepcounter{seqcount}
+ % the output
+ \edef\seqformat{\theseqcount}
+ % add subsection number
+ \ifsubsectionsineq
+ \let\tmp\seqformat
+ \edef\seqformat{\thesubsectioncount.\tmp}
+ \fi
+ % add section number
+ \ifsectionsineq
+ \let\tmp\seqformat
+ \edef\seqformat{\sectionprefix\thesectioncount.\tmp}
+ \fi
+ % define tag (for \label)
+ \xdef\tag{\seqformat}
+ % write number
+ \marginnote{\hfill(\seqformat)}
+}
+%% equation environment compatibility
+\def\equation{\hrefanchor$$\seqcount}
+\def\endequation{$$\@ignoretrue}
+
+
+%% figures
+% counter
+\newcounter{figcount}
+% booleans (write section or subsection in equation number)
+\newif\ifsectionsinfig
+\newif\ifsubsectionsinfig
+% width of figures
+\newlength\figwidth
+\setlength\figwidth\textwidth
+\addtolength\figwidth{-2.5cm}
+% caption
+\def\defcaption{
+ \long\def\caption##1{
+ \stepcounter{figcount}
+
+ % hyperref anchor
+ \hrefanchor
+
+ % the number of the figure
+ \edef\figformat{\thefigcount}
+ % add subsection number
+ \ifsubsectionsinfig
+ \let\tmp\figformat
+ \edef\figformat{\thesubsectioncount.\tmp}
+ \fi
+ % add section number
+ \ifsectionsinfig
+ \let\tmp\figformat
+ \edef\figformat{\sectionprefix\thesectioncount.\tmp}
+ \fi
+
+ % define tag (for \label)
+ \xdef\tag{\figformat}
+
+ % write
+ \hfil fig \figformat: \parbox[t]{\figwidth}{\leavevmode\small##1}
+
+ % space
+ \par\bigskip
+ }
+}
+%% short caption: centered
+\def\captionshort#1{
+ \stepcounter{figcount}
+
+ % hyperref anchor
+ \hrefanchor
+
+ % the number of the figure
+ \edef\figformat{\thefigcount}
+ % add section number
+ \ifsectionsinfig
+ \let\tmp\figformat
+ \edef\figformat{\sectionprefix\thesectioncount.\tmp}
+ \fi
+
+ % define tag (for \label)
+ \xdef\tag{\figformat}
+
+ % write
+ \hfil fig \figformat: {\small#1}
+
+ %space
+ \par\bigskip
+}
+
+%% environment
+\def\figure{
+ \par
+ \vfil\penalty100\vfilneg
+ \bigskip
+}
+\def\endfigure{
+ \par
+ \vfil\penalty100\vfilneg
+ \bigskip
+}
+
+
+%% start appendices
+\def\appendix{
+ \vfill
+ \pagebreak
+
+ % counter
+ \setcounter{sectioncount}0
+
+ % prefix
+ \def\sectionprefix{A}
+
+ % write
+ {\bf \LARGE Appendices}\par\penalty10000\bigskip\penalty10000
+
+ % add a mention in the table of contents
+ \iftoc
+ \immediate\write\tocoutput{\noexpand\tocappendices}\penalty10000
+ \fi
+
+ %% uncomment for new page for each appendix
+ %\def\seqskip{\vfill\pagebreak}
+}
+
+
+%% bibliography
+% size of header
+\newlength\bibheader
+\def\thebibliography#1{
+ \hrefanchor
+
+ % add a mention in the table of contents
+ \iftoc
+ % save lncount in aux variable which is written to toc
+ \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@references\endcsname{\thelncount}}
+ \write\tocoutput{\noexpand\tocreferences{\thepage}}\penalty10000
+ \fi
+
+ % write
+ {\bf \LARGE References}\par\penalty10000\bigskip\penalty10000
+ % width of header
+ \settowidth\bibheader{[#1]}
+ \leftskip\bibheader
+}
+% end environment
+\def\endthebibliography{
+ \par\leftskip0pt
+}
+
+%% bibitem command
+\def\bibitem[#1]#2{%
+ \hrefanchor%
+ \outdef{label@cite#2}{#1}%
+ \hskip-\bibheader%
+ \makebox[\bibheader]{\cite{#2}\hfill}%
+}
+
+%% cite command (adapted from latex.ltx)
+% @tempswa is set to true when there is an optional argument
+\newif\@tempswa
+\def\cite{%
+ % check whether there is an optional argument (if there is none, add on empty '[]')
+ \@ifnextchar [{\@tempswatrue\@citex}{\@tempswafalse\@citex[]}%
+}
+% command with optional argument
+\def\@citex[#1]#2{\leavevmode%
+ % initialize loop
+ \let\@citea\@empty%
+ % format
+ \@cite{%
+ % loop over ',' separated list
+ \@for\@citeb:=#2\do{%
+ % text to add at each iteration of the loop (separator between citations)
+ \@citea\def\@citea{,\ }%
+ % add entry to citelist
+ \@writecitation{\@citeb}%
+ \ref{cite\@citeb}%
+ }%
+ }%
+ % add optional argument text (as an argument to '\@cite')
+ {#1}%
+}
+\def\@cite#1#2{%
+ [#1\if@tempswa , #2\fi]%
+}
+%% add entry to citelist after checking it has not already been added
+\def\@writecitation#1{%
+ \ifcsname if#1cited\endcsname%
+ \else%
+ \expandafter\newif\csname if#1cited\endcsname%
+ \immediate\write\@auxout{\string\citation{#1}}%
+ \fi%
+}
+
+%% table of contents
+% boolean
+\newif\iftoc
+\def\tableofcontents{
+ {\bf \large Table of contents:}\par\penalty10000\bigskip\penalty10000
+
+ % copy content from file
+ \IfFileExists{\jobname.toc}{\input{\jobname.toc}}{{\tt error: table of contents missing}}
+
+ % open new toc
+ \newwrite\tocoutput
+ \immediate\openout\tocoutput=\jobname.toc
+
+ \toctrue
+}
+%% close file
+\AtEndDocument{
+ % close toc
+ \iftoc
+ \immediate\closeout\tocoutput
+ \fi
+}
+
+
+%% fill line with dots
+\def\leaderfill{\leaders\hbox to 1em {\hss. \hss}\hfill}
+
+%% same as sectionprefix
+\def\tocsectionprefix{}
+
+%% toc formats
+\newcounter{tocsectioncount}
+\def\tocsection #1#2{
+ \stepcounter{tocsectioncount}
+ \setcounter{tocsubsectioncount}{0}
+ \setcounter{tocsubsubsectioncount}{0}
+ % write
+ \smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par
+}
+\newcounter{tocsubsectioncount}
+\def\tocsubsection #1#2{
+ \stepcounter{tocsubsectioncount}
+ \setcounter{tocsubsubsectioncount}{0}
+ % write
+ {\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
+\newcounter{tocsubsubsectioncount}
+\def\tocsubsubsection #1#2{
+ \stepcounter{tocsubsubsectioncount}
+ % write
+ {\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
+\def\tocappendices{
+ \medskip
+ \setcounter{tocsectioncount}0
+ {\bf Appendices}\par
+ \smallskip
+ \def\tocsectionprefix{A}
+}
+\def\tocreferences#1{
+ \medskip
+ {\hyperlink{ln.\csname toc@references\endcsname}{{\color{blue}\bf References}\leaderfill#1}}\par
+ \smallskip
+}
+
+
+%% definitions that must be loaded at begin document
+\let\ian@olddocument\document
+\def\document{
+ \ian@olddocument
+
+ \deflabel
+ \defcaption
+}
+
+%% end
+\ian@defaultoptions
+\endinput
diff --git a/docs/nstrophy_doc/libs/iantheo.sty b/docs/nstrophy_doc/libs/iantheo.sty
new file mode 100644
index 0000000..d33a93d
--- /dev/null
+++ b/docs/nstrophy_doc/libs/iantheo.sty
@@ -0,0 +1,162 @@
+%%
+%% iantheorem package:
+%% Ian's customized theorem command
+%%
+
+%% boolean to signal that this package was loaded
+\newif\ifiantheo
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{iantheo}[2016/11/10]
+
+%% options
+\newif\ifsectionintheo
+\DeclareOption{section_in_theo}{\sectionintheotrue}
+\DeclareOption{no_section_in_theo}{\sectionintheofalse}
+\newif\ifsubsectionintheo
+\DeclareOption{subsection_in_theo}{\subsectionintheotrue}
+\DeclareOption{no_subsection_in_theo}{\subsectionintheofalse}
+
+\def\iantheo@defaultoptions{
+ \ExecuteOptions{section_in_theo, no_subsection_in_theo}
+ \ProcessOptions
+
+ %%% reset at every new section
+ \ifsectionintheo
+ \let\iantheo@oldsection\section
+ \gdef\section{\setcounter{theocount}{0}\iantheo@oldsection}
+ \fi
+
+ %% reset at every new subsection
+ \ifsubsectionintheo
+ \let\iantheo@oldsubsection\subsection
+ \gdef\subsection{\setcounter{theocount}{0}\iantheo@oldsubsection}
+ \fi
+}
+
+
+%% delimiters
+\def\delimtitle#1{
+ \par%
+ \leavevmode%
+ \raise.3em\hbox to\hsize{%
+ \lower0.3em\hbox{\vrule height0.3em}%
+ \hrulefill%
+ \ \lower.3em\hbox{#1}\ %
+ \hrulefill%
+ \lower0.3em\hbox{\vrule height0.3em}%
+ }%
+ \par\penalty10000%
+}
+
+%% callable by ref
+\def\delimtitleref#1{
+ \par%
+%
+ \ifdefined\ianclass%
+ % hyperref anchor%
+ \hrefanchor%
+ % define tag (for \label)%
+ \xdef\tag{#1}%
+ \fi%
+%
+ \leavevmode%
+ \raise.3em\hbox to\hsize{%
+ \lower0.3em\hbox{\vrule height0.3em}%
+ \hrulefill%
+ \ \lower.3em\hbox{\bf #1}\ %
+ \hrulefill%
+ \lower0.3em\hbox{\vrule height0.3em}%
+ }%
+ \par\penalty10000%
+}
+
+%% no title
+\def\delim{
+ \par%
+ \leavevmode\raise.3em\hbox to\hsize{%
+ \lower0.3em\hbox{\vrule height0.3em}%
+ \hrulefill%
+ \lower0.3em\hbox{\vrule height0.3em}%
+ }%
+ \par\penalty10000%
+}
+
+%% end delim
+\def\enddelim{
+ \par\penalty10000%
+ \leavevmode%
+ \raise.3em\hbox to\hsize{%
+ \vrule height0.3em\hrulefill\vrule height0.3em%
+ }%
+ \par%
+}
+
+
+%% theorem
+% counter
+\newcounter{theocount}
+% booleans (write section or subsection in equation number)
+\def\theo#1{
+ \stepcounter{theocount}
+ \ifdefined\ianclass
+ % hyperref anchor
+ \hrefanchor
+ \fi
+ % the number
+ \def\formattheo{\thetheocount}
+ % add subsection number
+ \ifsubsectionintheo
+ \let\tmp\formattheo
+ \edef\formattheo{\thesubsectioncount.\tmp}
+ \fi
+ % add section number
+ \ifsectionintheo
+ \let\tmp\formattheo
+ \edef\formattheo{\sectionprefix\thesectioncount.\tmp}
+ \fi
+ % define tag (for \label)
+ \xdef\tag{\formattheo}
+ % write
+ \delimtitle{\bf #1 \formattheo}
+}
+\let\endtheo\enddelim
+%% theorem headers with name
+\def\theoname#1#2{
+ \theo{#1}\hfil({\it #2})\par\penalty10000\medskip%
+}
+
+
+%% qed symbol
+\def\qedsymbol{$\square$}
+\def\qed{\penalty10000\hfill\penalty10000\qedsymbol}
+
+
+%% compatibility with article class
+\ifdefined\ianclasstrue
+ \relax
+\else
+ \def\thesectioncount{\thesection}
+ \def\thesubsectioncount{\thesubsection}
+ \def\sectionprefix{}
+\fi
+
+
+%% prevent page breaks after displayed equations
+\newcount\prevpostdisplaypenalty
+\def\nopagebreakaftereq{
+ \prevpostdisplaypenalty=\postdisplaypenalty
+ \postdisplaypenalty=10000
+}
+%% back to previous value
+\def\restorepagebreakaftereq{
+ \postdisplaypenalty=\prevpostdisplaypenalty
+}
+
+
+%% end
+\iantheo@defaultoptions
+\endinput
diff --git a/docs/nstrophy_doc/libs/largearray.sty b/docs/nstrophy_doc/libs/largearray.sty
new file mode 100644
index 0000000..ad5753b
--- /dev/null
+++ b/docs/nstrophy_doc/libs/largearray.sty
@@ -0,0 +1,19 @@
+%%
+%% largearray package:
+%% Array spanning the entire line
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{largearray}[2016/11/10]
+
+\RequirePackage{array}
+
+%% array spanning the entire line
+\newlength\largearray@width
+\setlength\largearray@width\textwidth
+\addtolength\largearray@width{-10pt}
+\def\largearray{\begin{array}{@{}>{\displaystyle}l@{}}\hphantom{\hspace{\largearray@width}}\\[-.5cm]}
+\def\endlargearray{\end{array}}
diff --git a/docs/nstrophy_doc/libs/point.sty b/docs/nstrophy_doc/libs/point.sty
new file mode 100644
index 0000000..9f983dc
--- /dev/null
+++ b/docs/nstrophy_doc/libs/point.sty
@@ -0,0 +1,106 @@
+%%
+%% Points package:
+%% \point commands
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{point}[2017/06/13]
+
+%% options
+\newif\ifresetatsection
+\DeclareOption{reset_at_section}{\resetatsectiontrue}
+\DeclareOption{no_reset_at_section}{\resetatsectionfalse}
+\newif\ifresetatsubsection
+\DeclareOption{reset_at_subsection}{\resetatsubsectiontrue}
+\DeclareOption{no_reset_at_subsection}{\resetatsubsectionfalse}
+\newif\ifresetattheo
+\DeclareOption{reset_at_theo}{\resetattheotrue}
+\DeclareOption{no_reset_at_theo}{\resetattheofalse}
+
+\def\point@defaultoptions{
+ \ExecuteOptions{reset_at_section, reset_at_subsection, no_reset_at_theo}
+ \ProcessOptions
+
+ %% reset at every new section
+ \ifresetatsection
+ \let\point@oldsection\section
+ \gdef\section{\resetpointcounter\point@oldsection}
+ \fi
+ %% reset at every new subsection
+ \ifresetatsubsection
+ \let\point@oldsubsection\subsection
+ \gdef\subsection{\resetpointcounter\point@oldsubsection}
+ \fi
+
+ %% reset at every new theorem
+ \ifresetattheo
+ \ifdefined\iantheotrue
+ \let\point@oldtheo\theo
+ \gdef\theo{\resetpointcounter\point@oldtheo}
+ \fi
+ \fi
+}
+
+
+%% point
+% counter
+\newcounter{pointcount}
+\def\point{
+ \stepcounter{pointcount}
+ \setcounter{subpointcount}{0}
+ % hyperref anchor (only if the class is 'ian')
+ \ifdefined\ifianclass
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\thepointcount}
+ \fi
+ % header
+ \indent{\bf \thepointcount\ - }
+}
+
+%% subpoint
+% counter
+\newcounter{subpointcount}
+\def\subpoint{
+ \stepcounter{subpointcount}
+ \setcounter{subsubpointcount}0
+ % hyperref anchor (only if the class is 'ian')
+ \ifdefined\ifianclass
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\thepointcount-\thesubpointcount}
+ \fi
+ % header
+ \indent\hskip.5cm{\bf \thepointcount-\thesubpointcount\ - }
+}
+
+%% subsubpoint
+% counter
+\newcounter{subsubpointcount}
+\def\subsubpoint{
+ \stepcounter{subsubpointcount}
+ % hyperref anchor (only if the class is 'ian')
+ \ifdefined\ifianclass
+ \hrefanchor
+ % define tag (for \label)
+ \xdef\tag{\thepointcount-\thesubpointcount-\thesubsubpointcount}
+ \fi
+ \indent\hskip1cm{\bf \thepointcount-\thesubpointcount-\thesubsubpointcount\ - }
+}
+
+
+%% reset point counters
+\def\resetpointcounter{
+ \setcounter{pointcount}{0}
+ \setcounter{subpointcount}{0}
+ \setcounter{subsubpointcount}{0}
+}
+
+
+
+%% end
+\point@defaultoptions
+\endinput
diff --git a/docs/nstrophy_doc/nstrophy_doc.tex b/docs/nstrophy_doc/nstrophy_doc.tex
new file mode 100644
index 0000000..fa5fff0
--- /dev/null
+++ b/docs/nstrophy_doc/nstrophy_doc.tex
@@ -0,0 +1,146 @@
+\documentclass{ian}
+
+\usepackage{largearray}
+
+\begin{document}
+
+\hbox{}
+\hfil{\bf\LARGE
+{\tt nstrophy}
+}
+\vfill
+
+\tableofcontents
+
+\vfill
+\eject
+
+\setcounter{page}1
+\pagestyle{plain}
+
+\section{Description of the computation}
+\subsection{Irreversible equation}
+\indent Consider the {\it irreversible} Navier-Stokes equation in 2 dimensions
+\begin{equation}
+ \partial_tu=\nu\Delta u+g-\nabla w-(u\cdot\nabla)u,\quad
+ \nabla\cdot u=0
+ \label{ins}
+\end{equation}
+in which $g$ is the forcing term and $w$ is the pressure.
+We take periodic boundary conditions, so, at every given time, $u(t,\cdot)$ is a function on the unit torus $\mathbb T^2:=\mathbb R^2/\mathbb Z^2$. We represent $u(t,\cdot)$ using its Fourier series
+\begin{equation}
+ \hat u_k(t):=\int_{\mathbb T^2}dx\ e^{2i\pi kx}u(t,x)
+\end{equation}
+for $k\in\mathbb Z^2$, and rewrite~\-(\ref{ins}) as
+\begin{equation}
+ \partial_t\hat u_k=
+ -4\pi^2\nu k^2\hat u_k+\hat g_k-2i\pi k\hat w_k
+ -2i\pi\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+ (q\cdot\hat u_p)\hat u_q
+ ,\quad
+ k\cdot\hat u_k=0
+ \label{ins_k}
+\end{equation}
+We then reduce the equation to a scalar one, by writing
+\begin{equation}
+ \hat u_k=\frac{2i\pi k^\perp}{|k|}\hat\varphi_k\equiv\frac{2i\pi}{|k|}(-k_y\hat\varphi_k,k_x\hat\varphi_k)
+\end{equation}
+in terms of which, multiplying both sides of the equation by $\frac{k^\perp}{|k|}$,
+\begin{equation}
+ \partial_t\hat \varphi_k=
+ -4\pi^2\nu k^2\hat \varphi_k+\frac{k^\perp}{2i\pi|k|}\cdot\hat g_k
+ -\frac{2i\pi}{|k|}\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+ (q\cdot p^\perp)(k^\perp\cdot q^\perp)\hat\varphi_p\hat\varphi_q
+ \label{ins_k}
+\end{equation}
+which, since $q\cdot p^\perp=q\cdot(k^\perp-q^\perp)=q\cdot k^\perp$, is
+\begin{equation}
+ \partial_t\hat \varphi_k=
+ -4\pi^2\nu k^2\hat \varphi_k+\frac{k^\perp}{2i\pi|k|}\cdot\hat g_k
+ +\frac{4\pi^2}{|k|}\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+ (q\cdot k^\perp)(k\cdot q)\frac{\hat\varphi_p\hat\varphi_q}{|p||q|}
+ .
+ \label{ins_k}
+\end{equation}
+We truncate the Fourier modes and assume that $\hat\varphi_k=0$ if $|k_1|>K_1$ or $|k_2|>K_2$. Let
+\begin{equation}
+ \mathcal K:=\{(k_1,k_2),\ |k_1|\leqslant K_1,\ |k_2|\leqslant K_2\}
+ .
+\end{equation}
+\bigskip
+
+\point{\bf FFT}. We compute the last term in~\-(\ref{ins_k})
+\begin{equation}
+ T(\hat\varphi,k):=
+ \sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+ \frac{\hat\varphi_p}{|p|}
+ (q\cdot k^\perp)(k\cdot q)\frac{\hat\varphi_q}{|q|}
+\end{equation}
+using a fast Fourier transform, defined as
+\begin{equation}
+ \mathcal F(f)(n):=\sum_{m\in\mathcal N}e^{-\frac{2i\pi}{N_1}m_1n_1-\frac{2i\pi}{N_2}m_2n_2}f(m_1,m_2)
+\end{equation}
+where
+\begin{equation}
+ \mathcal N:=\{(n_1,n_2),\ 0\leqslant n_1\leqslant N_1,\ 0\leqslant n_2\leqslant N_2\}
+\end{equation}
+for some fixed $N_1,N_2$. The transform is inverted by
+\begin{equation}
+ \frac1{N_1N_2}\mathcal F^*(\mathcal F(f))(n)=f(n)
+\end{equation}
+in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase.
+\bigskip
+
+\indent The condition $p+q=k$ can be rewritten as
+\begin{equation}
+ T(\hat\varphi,k)
+ =
+ \sum_{p,q\in\mathcal K}
+ \frac1{N_1N_2}
+ \sum_{n\in\mathcal N}e^{-\frac{2i\pi}{N_1}n_1(p_1+q_1-k_1)-\frac{2i\pi}{N_2}n_2(p_2+q_2-k_2)}
+ \frac{\hat\varphi_p}{|p|}
+ (q\cdot k^\perp)(k\cdot q)\frac{\hat\varphi_q}{|q|}
+\end{equation}
+provided
+\begin{equation}
+ N_i>4K_i.
+\end{equation}
+Indeed, $\sum_{n_i=0}^{N_i}e^{-\frac{2i\pi}{N_i}n_im_i}$ vanishes unless $m_i=0\%N_i$ (in which $\%N_i$ means `modulo $N_i$'), and, if $p,q\in\mathcal K$, then $|p_i+q_i|\leqslant2K_i$. Therefore,
+\begin{equation}
+ T(\hat\varphi,k)
+ =
+ \frac1{N_1N_2}
+ \mathcal F^*\left(
+ \mathcal F(|p|^{-1}\hat\varphi_p)(n)
+ \mathcal F((|q|^{-1}(q\cdot k^\perp)(k\cdot q)\hat\varphi_q)_q)(n)
+ \right)(k)
+\end{equation}
+which we expand
+\begin{equation}
+ T(\hat\varphi,k)
+ =
+ \textstyle
+ \frac{k_x^2-k_y^2}{N_1N_2}
+ \mathcal F^*\left(
+ \mathcal F\left(\frac{\hat\varphi_p}{|p|}\right)(n)
+ \mathcal F\left(\frac{q_xq_y}{|q|}\hat\varphi_q\right)(n)
+ \right)(k)
+ -
+ \frac{k_xk_y}{N_1N_2}
+ \mathcal F^*\left(
+ \mathcal F\left(\frac{\hat\varphi_p}{|p|}\right)(n)
+ \mathcal F\left(\frac{q_x^2-q_y^2}{|q|}\hat\varphi_q\right)(n)
+ \right)(k)
+\end{equation}
+
+
+\vfill
+\eject
+
+\begin{thebibliography}{WWW99}
+\small
+\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}
+\end{thebibliography}
+
+
+\end{document}
diff --git a/src/main.c b/src/main.c
new file mode 100644
index 0000000..9420e90
--- /dev/null
+++ b/src/main.c
@@ -0,0 +1,282 @@
+#define VERSION "0.0"
+
+#include <math.h>
+#include <complex.h>
+#include <fftw3.h>
+#include <string.h>
+#include <stdlib.h>
+#include "navier-stokes.h"
+
+// usage message
+int print_usage();
+// read command line arguments
+int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nsteps, unsigned int* computation_nr);
+
+// compute enstrophy as a function of time in the I-NS equation
+int enstrophy(ns_params params, unsigned int Nsteps);
+
+
+#define COMPUTATION_ENSTROPHY 1
+int main (int argc, const char* argv[]){
+ ns_params params;
+ unsigned int nsteps;
+ int ret;
+ unsigned int computation_nr;
+
+ // default computation: phase diagram
+ computation_nr=COMPUTATION_ENSTROPHY;
+
+ // read command line arguments
+ ret=read_args(argc, argv, &params, &nsteps, &computation_nr);
+ if(ret<0){
+ return(-1);
+ }
+ if(ret>0){
+ return(0);
+ }
+
+ // enstrophy
+ if(computation_nr==COMPUTATION_ENSTROPHY){
+ enstrophy(params, nsteps);
+ }
+
+ return(0);
+}
+
+// usage message
+int print_usage(){
+ fprintf(stderr, "usage:\n nstrophy enstrophy [-h timestep] [-K modes] [-v] [-N nsteps]\n\n nstrophy -V [-v]\n\n");
+ return(0);
+}
+
+// read command line arguments
+#define CP_FLAG_TIMESTEP 1
+#define CP_FLAG_NSTEPS 2
+#define CP_FLAG_MODES 2
+#define CP_FLAG_NU 3
+int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nsteps, unsigned int* computation_nr){
+ int i;
+ int ret;
+ // temporary int
+ int tmp_int;
+ // temporary unsigned int
+ unsigned int tmp_uint;
+ // temporary double
+ double tmp_double;
+ // pointers
+ char* ptr;
+ // flag that indicates what argument is being read
+ int flag=0;
+ // print version and exit
+ char Vflag=0;
+
+ // defaults
+ /*
+ params->h=6.103515625e-05;
+ params->K=16;
+ *nsteps=16777216;
+ params->nu=4.9632717887631524e-05;
+ */
+ params->h=1e-5;
+ params->K=16;
+ *nsteps=10000000;
+ params->nu=1e-4;
+
+ // loop over arguments
+ for(i=1;i<argc;i++){
+ // flag
+ if(argv[i][0]=='-'){
+ for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
+ switch(*ptr){
+ // timestep
+ case 'h':
+ flag=CP_FLAG_TIMESTEP;
+ break;
+ // nsteps
+ case 'N':
+ flag=CP_FLAG_NSTEPS;
+ break;
+ // modes
+ case 'K':
+ flag=CP_FLAG_MODES;
+ break;
+ // friction
+ case 'n':
+ flag=CP_FLAG_NU;
+ break;
+ // print version
+ case 'V':
+ Vflag=1;
+ break;
+ default:
+ fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
+ print_usage();
+ return(-1);
+ break;
+ }
+ }
+ }
+ // timestep
+ else if(flag==CP_FLAG_TIMESTEP){
+ ret=sscanf(argv[i],"%lf",&tmp_double);
+ if(ret!=1){
+ fprintf(stderr, "error: '-h' should be followed by a double\n got '%s'\n",argv[i]);
+ return(-1);
+ }
+ params->h=tmp_double;
+ flag=0;
+ }
+ // nsteps
+ else if(flag==CP_FLAG_NSTEPS){
+ ret=sscanf(argv[i],"%u",&tmp_uint);
+ if(ret!=1){
+ fprintf(stderr, "error: '-N' should be followed by an unsigned int\n got '%s'\n",argv[i]);
+ return(-1);
+ }
+ *nsteps=tmp_uint;
+ flag=0;
+ }
+ // nsteps
+ else if(flag==CP_FLAG_MODES){
+ ret=sscanf(argv[i],"%d",&tmp_int);
+ if(ret!=1){
+ fprintf(stderr, "error: '-K' should be followed by an int\n got '%s'\n",argv[i]);
+ return(-1);
+ }
+ params->K=tmp_uint;
+ flag=0;
+ }
+ // friction
+ else if(flag==CP_FLAG_TIMESTEP){
+ ret=sscanf(argv[i],"%lf",&tmp_double);
+ if(ret!=1){
+ fprintf(stderr, "error: '-n' should be followed by a double\n got '%s'\n",argv[i]);
+ return(-1);
+ }
+ params->nu=tmp_double;
+ flag=0;
+ }
+ // computation to run
+ else{
+ if(strcmp(argv[i], "enstrophy")==0){
+ *computation_nr=COMPUTATION_ENSTROPHY;
+ }
+ else{
+ fprintf(stderr, "error: unrecognized computation: '%s'\n",argv[i]);
+ print_usage();
+ return(-1);
+ }
+ flag=0;
+ }
+ }
+
+ // print version and exit
+ if(Vflag==1){
+ printf("nstrophy " VERSION "\n");
+ return(1);
+ }
+
+ return(0);
+}
+
+// compute enstrophy as a function of time in the I-NS equation
+int enstrophy(ns_params params, unsigned int Nsteps){
+ _Complex double* u;
+ _Complex double* tmp1;
+ _Complex double* tmp2;
+ _Complex double* tmp3;
+ _Complex double alpha;
+ _Complex double avg;
+ unsigned int t;
+ int kx,ky;
+ fft_vects fft_vects;
+
+ // sizes
+ params.S=2*params.K+1;
+ params.N=4*params.K+1;
+
+ // velocity field
+ u=calloc(sizeof(_Complex double),params.S*params.S);
+ params.g=calloc(sizeof(_Complex double),params.S*params.S);
+ // allocate tmp vectors for computation
+ tmp1=calloc(sizeof(_Complex double),params.S*params.S);
+ tmp2=calloc(sizeof(_Complex double),params.S*params.S);
+ tmp3=calloc(sizeof(_Complex double),params.S*params.S);
+
+ // initial value
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ //u[KLOOKUP(kx,ky,params.S)]=kx*kx*ky*ky*exp(-(kx*kx+ky*ky));
+ if((kx==1 && ky==0) || (kx==-1 && ky==0)){
+ u[KLOOKUP(kx,ky,params.S)]=1;
+ }
+ else{
+ u[KLOOKUP(kx,ky,params.S)]=0;
+ }
+ }
+ }
+
+ // driving force
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ //params.g[KLOOKUP(kx,ky,params.S)]=sqrt(kx*kx*ky*ky)*exp(-(kx*kx+ky*ky));
+ if((kx==2 && ky==-1) || (kx==-2 && ky==1)){
+ params.g[KLOOKUP(kx,ky,params.S)]=1;
+ }
+ else{
+ params.g[KLOOKUP(kx,ky,params.S)]=0;
+ }
+ }
+ }
+
+
+ // prepare vectors for fft
+ fft_vects.fft1=fftw_malloc(sizeof(fftw_complex)*params.N*params.N);
+ fft_vects.fft1_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.fft1, fft_vects.fft1, FFTW_FORWARD, FFTW_MEASURE);
+ fft_vects.fft2=fftw_malloc(sizeof(fftw_complex)*params.N*params.N);
+ fft_vects.fft2_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.fft2, fft_vects.fft2, FFTW_FORWARD, FFTW_MEASURE);
+ fft_vects.invfft=fftw_malloc(sizeof(fftw_complex)*params.N*params.N);
+ fft_vects.invfft_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.invfft, fft_vects.invfft, FFTW_BACKWARD, FFTW_MEASURE);
+
+ // init running average
+ avg=0;
+
+ // iterate
+ for(t=0;t<Nsteps;t++){
+ ins_step(u, params, fft_vects, tmp1, tmp2, tmp3);
+ alpha=compute_alpha(u, params);
+
+ // to avoid errors building up in imaginary part
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ u[KLOOKUP(kx,ky,params.S)]=__real__ u[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+
+ // running average
+ if(t>0){
+ avg=avg-(avg-alpha)/t;
+ }
+
+ if(t>0 && t%1000==0){
+ fprintf(stderr,"%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha);
+ printf("%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha);
+ }
+ }
+
+ // free memory
+ fftw_destroy_plan(fft_vects.fft1_plan);
+ fftw_destroy_plan(fft_vects.fft2_plan);
+ fftw_destroy_plan(fft_vects.invfft_plan);
+ fftw_free(fft_vects.fft1);
+ fftw_free(fft_vects.fft2);
+ fftw_free(fft_vects.invfft);
+
+ free(tmp3);
+ free(tmp2);
+ free(tmp1);
+ free(params.g);
+ free(u);
+
+ return(0);
+}
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
new file mode 100644
index 0000000..3def963
--- /dev/null
+++ b/src/navier-stokes.c
@@ -0,0 +1,299 @@
+#include "navier-stokes.h"
+#include <math.h>
+
+#define M_PI 3.14159265358979323846
+
+#define CHK 1
+
+// next time step for Irreversible Navier-Stokes equation
+int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3){
+ int kx,ky;
+
+ // k1
+ ins_rhs(tmp1, u, params, vects);
+ // add to output
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ tmp3[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/6*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+
+ // u+h*k1/2
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/2*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+ // k2
+ ins_rhs(tmp1, tmp2, params, vects);
+ // add to output
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ tmp3[KLOOKUP(kx,ky,params.S)]+=params.h/3*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+
+ // u+h*k2/2
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/2*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+ // k3
+ ins_rhs(tmp1, tmp2, params, vects);
+ // add to output
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ tmp3[KLOOKUP(kx,ky,params.S)]+=params.h/3*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+
+ // u+h*k3
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+ // k4
+ ins_rhs(tmp1, tmp2, params, vects);
+ // add to output
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ u[KLOOKUP(kx,ky,params.S)]=tmp3[KLOOKUP(kx,ky,params.S)]+params.h/6*tmp1[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+
+ return(0);
+}
+
+// right side of Irreversible Navier-Stokes equation
+int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects){
+ int kx,ky;
+
+#if CHK==0
+ // F(u/|p|)*F(q1*q2*u/|q|)
+ // init to 0
+ for(kx=0; kx<params.N*params.N; kx++){
+ vects.fft1[kx]=0;
+ vects.fft2[kx]=0;
+ vects.invfft[kx]=0;
+ }
+ // fill modes
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ if(kx!=0 || ky!=0){
+ vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+ vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*ky*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+ }
+ }
+ }
+ // fft
+ fftw_execute(vects.fft1_plan);
+ fftw_execute(vects.fft2_plan);
+ // write to invfft
+ for(kx=-2*params.K;kx<=2*params.K;kx++){
+ for(ky=-2*params.K;ky<=2*params.K;ky++){
+ vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+ }
+ }
+ // inverse fft
+ fftw_execute(vects.invfft_plan);
+
+ // write out
+ for(kx=0; kx<params.S*params.S; kx++){
+ out[kx]=0;
+ }
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ if(kx!=0 || ky!=0){
+ out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*kx-ky*ky)/params.N/params.N;
+ }
+ }
+ }
+
+ // F(u/|p|)*F((q1*q1-q2*q2)*u/|q|)
+ // init to 0
+ for(kx=0; kx<params.N*params.N; kx++){
+ vects.fft2[kx]=0;
+ vects.invfft[kx]=0;
+ }
+ // fill modes
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ if(kx!=0 || ky!=0){
+ vects.fft2[KLOOKUP(kx,ky,params.N)]=(kx*kx-ky*ky)*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+ }
+ }
+ }
+ // fft
+ fftw_execute(vects.fft2_plan);
+ // write to invfft
+ for(kx=-2*params.K;kx<=2*params.K;kx++){
+ for(ky=-2*params.K;ky<=2*params.K;ky++){
+ vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+ }
+ }
+ // inverse fft
+ fftw_execute(vects.invfft_plan);
+
+ // write out
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ if(kx!=0 || ky!=0){
+ out[KLOOKUP(kx,ky,params.S)]=out[KLOOKUP(kx,ky,params.S)]-4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*ky)/params.N/params.N;
+ }
+ }
+ }
+#elif CHK == 1
+ // F(-p2/|p|*u)*F(q1*|q|*u)
+ // init to 0
+ for(kx=0; kx<params.N*params.N; kx++){
+ vects.fft1[kx]=0;
+ vects.fft2[kx]=0;
+ vects.invfft[kx]=0;
+ }
+ // fill modes
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ if(kx!=0 || ky!=0){
+ vects.fft1[KLOOKUP(kx,ky,params.N)]=-kx/sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
+ vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+ }
+ // fft
+ fftw_execute(vects.fft1_plan);
+ fftw_execute(vects.fft2_plan);
+ // write to invfft
+ for(kx=-2*params.K;kx<=2*params.K;kx++){
+ for(ky=-2*params.K;ky<=2*params.K;ky++){
+ vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)] - vects.fft1[KLOOKUP(ky,kx,params.N)]*vects.fft2[KLOOKUP(ky,kx,params.N)];
+ }
+ }
+
+ // inverse fft
+ fftw_execute(vects.invfft_plan);
+
+
+ // write out
+ for(kx=0; kx<params.S*params.S; kx++){
+ out[kx]=0;
+ }
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ if(kx!=0 || ky!=0){
+ out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu/params.S*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+2*M_PI/sqrt(kx*kx+ky*ky)/params.S*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
+ }
+ }
+ }
+
+#elif CHK==2
+ // F(u)*F(q1*u)
+ // init to 0
+ for(kx=0; kx<params.N*params.N; kx++){
+ vects.fft1[kx]=0;
+ vects.fft2[kx]=0;
+ vects.invfft[kx]=0;
+ }
+ // fill modes
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ // u_k
+ vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)];
+ // 2i\pi k_x u_k
+ vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*kx*u[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+ // fft
+ fftw_execute(vects.fft1_plan);
+ fftw_execute(vects.fft2_plan);
+ // write to invfft
+ for(kx=-2*params.K;kx<=2*params.K;kx++){
+ for(ky=-2*params.K;ky<=2*params.K;ky++){
+ vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+ }
+ }
+
+ // F(p1/p2*u)*F(q2*u)
+ // init to 0
+ for(kx=0; kx<params.N*params.N; kx++){
+ vects.fft1[kx]=0;
+ vects.fft2[kx]=0;
+ }
+ // fill modes
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ // k_x/k_y u_k
+ if(ky!=0){
+ vects.fft1[KLOOKUP(kx,ky,params.N)]=kx/ky*u[KLOOKUP(kx,ky,params.S)];
+ }
+ // 2i\pi k_y u_k
+ vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*ky*u[KLOOKUP(kx,ky,params.S)];
+ }
+ }
+ // fft
+ fftw_execute(vects.fft1_plan);
+ fftw_execute(vects.fft2_plan);
+ // write to invfft
+ for(kx=-2*params.K;kx<=2*params.K;kx++){
+ for(ky=-2*params.K;ky<=2*params.K;ky++){
+ vects.invfft[KLOOKUP(kx,ky,params.N)]+=-vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+ }
+ }
+
+ // inverse fft
+ fftw_execute(vects.invfft_plan);
+
+ /*
+ // check: convolution instead of fft
+ for(kx=0; kx<params.S*params.S; kx++){
+ out[kx]=0;
+ }
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ for(px=-params.K;px<=params.K;px++){
+ for(py=-params.K;py<=params.K;py++){
+ if(kx-px<=params.K && kx-px>=-params.K && ky-py<=params.K && ky-py>=-params.K){
+ out[KLOOKUP(kx,ky,params.S)]+=2*M_PI*I*(u[KLOOKUP(px,py,params.S)]*(kx-px)*u[KLOOKUP(kx-px,ky-py,params.S)]-(py==0?0:px/py*u[KLOOKUP(px,py,params.S)]*(ky-py)*u[KLOOKUP(kx-px,ky-py,params.S)]));
+ }
+ }
+ }
+ dd=(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])*(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])+(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)])*(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)]);
+ if(dd>1e-25){
+ printf("%d %d % .8e\n",kx,ky, dd);
+ }
+ }
+ }
+ */
+
+
+ // write out
+ for(kx=0; kx<params.S*params.S; kx++){
+ out[kx]=0;
+ }
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
+ }
+ }
+
+#endif
+ return(0);
+}
+
+
+// compute alpha
+_Complex double compute_alpha(_Complex double* u, ns_params params){
+ _Complex double num=0;
+ _Complex double denom=0;
+ int kx,ky;
+
+ for(kx=-params.K;kx<=params.K;kx++){
+ for(ky=-params.K;ky<=params.K;ky++){
+ denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]*conj(u[KLOOKUP(kx,ky,params.S)])*(1+(ky!=0?kx*kx/ky/ky:0));
+ num+=(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]*conj(params.g[KLOOKUP(kx,ky,params.S)])*(1+(ky!=0?kx*kx/ky/ky:0));
+ }
+ }
+
+ return(num/denom);
+}
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
new file mode 100644
index 0000000..cd093d7
--- /dev/null
+++ b/src/navier-stokes.h
@@ -0,0 +1,48 @@
+#ifndef NAVIERSTOKES_H
+#define NAVIERSTOKES_H
+
+#include <complex.h>
+#include <fftw3.h>
+
+// to extract elements from array of size S representing a function of momentum, use
+// array[KEXTRACT(kx,ky,size)]
+#define KLOOKUP(X,Y,S) (X>=0?X:S+X)*S+(Y>=0?Y:S+Y)
+
+
+// parameters for the NS equation
+typedef struct ns_params {
+ // number of modes
+ int K;
+ // 2*K+1
+ int S;
+ // 4*K+1
+ int N;
+ // forcing term
+ _Complex double* g;
+ // time step
+ double h;
+ // friction
+ double nu;
+} ns_params;
+
+// arrays on which the ffts are performed
+typedef struct fft_vects {
+ fftw_complex* fft1;
+ fftw_complex* fft2;
+ fftw_complex* invfft;
+ fftw_plan fft1_plan;
+ fftw_plan fft2_plan;
+ fftw_plan invfft_plan;
+} fft_vects;
+
+// next time step for Irreversible Navier-Stokes equation
+int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3);
+
+// right side of Irreversible Navier-Stokes equation
+int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects);
+
+// compute alpha
+_Complex double compute_alpha(_Complex double* u, ns_params params);
+
+#endif
+
generated by cgit v1.2.3-54-g00ecf (git 2.45.2) at 2024-07-03 13:20:16 +0000