From 01f47ace6756c28deb9ea0daaee3904ffa5ce9e0 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Thu, 11 Jan 2018 22:48:14 +0000 Subject: Initial commit --- Makefile | 54 +++ docs/nstrophy_doc/Makefile | 70 ++++ docs/nstrophy_doc/bibliography/conf.BBlog | 5 + docs/nstrophy_doc/libs/ian.cls | 659 ++++++++++++++++++++++++++++++ docs/nstrophy_doc/libs/iantheo.sty | 162 ++++++++ docs/nstrophy_doc/libs/largearray.sty | 19 + docs/nstrophy_doc/libs/point.sty | 106 +++++ docs/nstrophy_doc/nstrophy_doc.tex | 146 +++++++ src/main.c | 282 +++++++++++++ src/navier-stokes.c | 299 ++++++++++++++ src/navier-stokes.h | 48 +++ 11 files changed, 1850 insertions(+) create mode 100644 Makefile create mode 100644 docs/nstrophy_doc/Makefile create mode 100644 docs/nstrophy_doc/bibliography/conf.BBlog create mode 100644 docs/nstrophy_doc/libs/ian.cls create mode 100644 docs/nstrophy_doc/libs/iantheo.sty create mode 100644 docs/nstrophy_doc/libs/largearray.sty create mode 100644 docs/nstrophy_doc/libs/point.sty create mode 100644 docs/nstrophy_doc/nstrophy_doc.tex create mode 100644 src/main.c create mode 100644 src/navier-stokes.c create mode 100644 src/navier-stokes.h 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 +#include +#include +#include +#include +#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, ¶ms, &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;ih=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;t0){ + 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 + +#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.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 +#include + +// 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 + -- cgit v1.2.3-54-g00ecf