User Manual for 2D Code Kaifeng Chen Universit¨ at Heidelberg

1

Francis Chen

User Manual

Contents 1

The Structure

3

2

Constants

3

3

Input Interface

5 5 5 5

3.1 3.2 3.3 4

FFTW Solver 4.1 4.2 4.3

5

Private Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Input Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . FFTW Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Private Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Output Interface 5.1 5.2 5.3

7 7 7 9 10

Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Output Flags . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Output Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

6

Possible Problems

14

7

Source Code & Documents

14

2 of 14

Francis Chen

1

User Manual

The Structure

The 2 dimensional code of Burgers equation consists of 8 files and external packages FFTW. These 8 files are: 1. const.h: constants in the code. 2. input.h, input.cpp: header file and .cpp file for the input interface. 3. output.h, output.cpp: header file and .cpp file for the output interface. 4. solver FFTW.h, solver FFTW.cpp: header file and .cpp file for the process of solving the Burgers. 5. main.cpp: .cpp file that executes the code. All these independent files are divided into 3 classes (C++ term): Class Input, Class Output, Class Solver FFTW1 . The typical work flow of a single run is like this: Reading the input file (a .txt file) ⇒ solving Burgers equation using certain algorithms ⇒ generating outputs. The following context will give a more detailed description of all these files. In this pdf, all the classes and all the data types are in bold font like Input.

2

Constants

This section mainly discusses the file const.h. All the constants that might be used in the code is shown here. Here we use ”# define”2 to set all the constants: • VISCOSITY: viscosity set in the simulation, before rescaling. • TIME N: total time steps. • TIME STEP: time step size ∆t. • PI: circumference ratio, mathematical constant. • INPUT PATH: input path in one’s computer, when the code reads the input file, this path will be added before the name of the file. • OUTPUT PATH: output path in one’s computer, when the code generates an output file, this path will be added before the name of the file. • ENERGY OUTPUT: numbers of time steps in between calculations and output of the energy. • GENERATE OUTPUT: numbers of time steps in between generating two outputs. • THREADS: number of threads used in the parallel computing. All the constants in const.h can be changes according to one’s needs except for PI. Here is one example of setting INPUT PATH and OUTPUT PATH: 1A 2 It

convention in C++: the name of a class should begin with a capital letter is a convention in C++ to use capital letters here.

3 of 14

Francis Chen

User Manual

Example 1: If one’s input file is in ”D:\Heidelberg\2-D\ Executives\input\”, then one should assign this string to INPUT PATH. When the code begin to read the input, the console will tell you to enter the file name, then one only has to type input 1.txt. Here is a picture showing the example:

Figure 1: Example 1 The input path is actually: ”D:\Heidelberg\2-D\ Executives\input\input 1.txt”. So is the same with OUTPUT PATH. Remember one thing: the exact input and output path is ”INPUT PATH+input file name” or ”OUTPUT PATH+output file name”. Here is another example of ENERGY OUTPUT and GENERATE OUTPUT:

Example 2: If one choose ENERGY OUTPUT to be 10000, then after every 10000 steps the code will calculate the energy of the system (energy after rescaling) and put the result into E.txt. For example, the first time calculate energy is at t = 0, then the next time will be t = 10000. Here t means the steps instead of physical time. Here are some useful formulas one should know: physical time = T IM E ST EP × steps rescaled physical time = physical time × rescaling f actor

(1) (2)

4 of 14

Francis Chen

3

User Manual

Input Interface

In the code, input is set up to be a class, so if one want to create an input, one should use the interface like this: Input* aNewInput = new Input() ”aNewInput” will be a specific object, then one only has to type the input name and the code will automatically finish the input process. The details of the input will appear in the console, please look back to figure 1. These details include the number of grids in two directions and the number of threads. The following discussions in this section are all in input.cpp and input.h files.

3.1

Private Variables

The private variables in input.h are: • int numOfXGrid: number of grid points along x direction. • int numOfYGrid: number of grid points along y direction. • double** xVelocity: storing v(x, y, 0) of every grid point. • double** yVelocity: storing w(x, y, 0) of every grid point. The two arrays xVelocity and yVelocity are real arrays with the dimension: numOf XGrid × numOf Y Grid.

3.2

Functions

There are totally 6 functions in input.h: • Input(): constructor of Class Input. • ∼Input(): destructor of Class Input. • int getXGridNum(): return the value of numOfXGrid. • int getYGridNum(): return the value of numOfYGrid. • double getXVelocity(int i,int j): return the value of xVelocity[i][j]; • double getYVelocity(int i,int j): return the value of yVelocity[i][j]; The implementation of the 6 functions can be viewed in input.cpp.

3.3

Input Format

There are no requirements for the name of the input. The type of the input should end up with .txt or .dat. The input data are always generated by Matlab and then the code reads the data from the hard disk. Here are the requirements for the date inside the file: • the first line in the file should be two integers indicating the number of the grid points in x axis and in y axis.

5 of 14

Francis Chen

User Manual

• the following data are made up of 4 columns, the first two columns are the coordinates of the grid point, and the last two are the v(x, y, 0) and w(x, y, 0) of that point. • by now all these coordinates and velocities are before rescaling. • all the data should be separated by ’\t’ or ’\n’.

Example 3: here is an example of a typical input file:

Figure 2: Example 3 One can refer to the example Input in github repository.

6 of 14

Francis Chen

4

User Manual

FFTW Solver

If one wants to create an object from Class solver FFTW, the code is like this: Solver FFTW* aNewTest = new Solver FFTW() This only appears in main.cpp where the code is executed.

4.1

FFTW Package

Two functions in package FFTW are used in this class: • fftw plan dft r2c 2d(int n 0,int n 1, double *in, fftw complex *out, unsigned flags): do the FFT and get the output fftw complex *out. • fftw plan dft c2r 2d(int n 0,int n 1, fftw complex *in, double *out, unsigned flags): do the IFFT and get the output double *out. For the r2c FFT, the last dimension of the output complex array is only half of the initial size, which means even the real arrays are of dimensions: Nx × Ny , the complex arrays in fourier space are of dimension: Nx × (Ny /2 + 1). Therefore, the last dimension set in the code is y. FFTW already has the data type for the complex numbers:fftw complex. And if one wants to do the transforms, one can generate plans and then execute the plans to perform the transforms. Since the flag FFTW PRESERVE INPUT is not applicable in 2D, so the input array will be destroyed after every execution. To avoid this problem, some temporary arrays are used as the inputs. More information and details can be found in the user manual of FFTW in its homepage.

4.2

Private Variables

In solver FFTW.h, there are totally 32 variables: • int numOfXGrid: number of grid points along x direction. • int numOfYGrid: number of grid points along y direction. • double initE: initial energy of the system. See equation (5). • double**v: v(x, y). • double**w: w(x, y). • double**temp Velocity: temporary array to store velocities when doing FFT and IFFT. • double**firstD u: temporary arrays to store first order derivatives. • double**secondD u: temporary arrays to store second order derivatives. • double**v x: array for ∂v/∂x. • double**v y: array for ∂v/∂y. • double**w x: array for ∂w/∂x.

7 of 14

Francis Chen

User Manual

• double**w y: array for ∂w/∂y. • double**v x x: array for ∂ 2 v/∂x2 . • double**v y y: array for ∂ 2 v/∂y 2 . • double**w x x: array for ∂ 2 w/∂x2 . • double**w y y: array for ∂ 2 w/∂y 2 . • double**externalFx: array for external force along x direction. • double**externalFy: array for external force along y direction. • double***Adams v: array for Adams-Bashforth method for v(x, y). • double***Adams w: array for Adams-Bashforth method for w(x, y). • fftw complex**V: fourier transform for v(x, y). • fftw complex**W: fourier transform for w(x, y). • fftw complex*temp U: temporary array when doing the FFT and IFFT. • fftw complex*firstD U: temporary arrays to store fourier transform of first order derivatives. • fftw complex*secondD U: temporary arrays to store fourier transform of second order derivatives. • fftw complex**temp: temporary array when transforming 2D array to 1D in initializing the plan. • fftw plan plan r2c: plan for FFT. • fftw plan plan c2r: plan for IFFT. • fftw plan plan firstD: plan for calculating first order derivatives. • fftw plan plan secondD: plan for calculating second order derivatives. • ofstream energy: store the energy with respect to time with the file name E.txt. The time and the energy in this file are rescaled. • ofstream readMe: store the parameters set in the simulation with the file name readM e.txt. The following picture is an example of a readMe file:

Example 4:

Figure 3: Example 4

8 of 14

Francis Chen

4.3

User Manual

Functions

In solver FFTW.h, there are totally 8 functions, details of these functions can be viewed in solver FFTW.cpp. • Solver FFTW(): constructor of Class Solver FFTW. • ∼Solver FFTW(): destructor of Class Solver FFTW. • void burgersSolver FFTW(): function solving Burgers equation. • void firstDerivative(): function calculating first order derivatives. • void secondDerivative(): function calculating second order derivatives. • void adamsMethod(int t): doing Adams-Bashforth method in time step t. • void samplingForce(): sampling external force. By now it is a dummy function. • double calculateE(): calculating energy (before rescaling). One only needs to complete void samplingForce() if one wants to add external force to Burgers equation.

9 of 14

Francis Chen

5

User Manual

Output Interface

In the code, the output is a class. If one want to generate an output, the code is like this: Output* aNewOutput = new Output(int int 1, int int 2, double** array 1, double** array 2, int time, double energy, output Type FLAG) Here are some explanations of all the variables: • int int 1: int 1=numOfXGrid. • int int 2: int 2=numOfYGrid. • int time: the time step when this output is generated. • double** array 1: a two dimensional real array related to v(x, y, time). See Output Flags. • double** array 2: a two dimensional real array related to w(x, y, time). See Output Flags. • double energy: the initial energy, used for rescaling. • output Type FLAG: output type. See Output Flags.

5.1

Functions

There are only two functions in Class Output: • Output(int int 1, int int 2, double** array 1, double** array 2, int time, double energy, output Type FLAG): the constructor of Class Output. • ∼Output(): the destructor of Class Output. Details of the functions can be viewed in output.cpp.

5.2

Output Flags

To make the output interface powerful, here I set five flags to choose. Each flag stands for a different type of output: • VELOCITY: generating an output for the velocities v, w at time step t. Output* aNewOutput = new Output(numOfXGrid, numOfYGrid, v, w, t, initE, VELOCITY) • DERIVATIVEv: generating an output for the derivatives v x, v y at time step t. Output* aNewOutput = new Output(numOfXGrid, numOfYGrid, v x, v y, t, initE, DERIVATIVEv) • DERIVATIVEw: generating an output for the derivatives w x, w y time step t. Output* aNewOutput = new Output(numOfXGrid, numOfYGrid, w x, w y, t, initE, DERIVATIVEw) • DDERIVATIVEv: generating an output for the derivatives v x x, v y y time step t. 10 of 14

Francis Chen

User Manual

Output* aNewOutput = new Output(numOfXGrid, numOfYGrid, v x x, v y y, t, initE, DDERIVATIVEv) • DDERIVATIVEw: generating an output for the derivatives w x x, w y y time step t. Output* aNewOutput = new Output(numOfXGrid, numOfYGrid, w x x, w y y, t, initE, DDERIVATIVEw)

Example 5: here is a example in the code to generate outputs:

Figure 4: Example 5 The rescaling of the system is as follows (see the summary): 1 Lx 1 y0 = y · Lx

x0 = x ·

1 v0 = v · p E · Ly /Lx 1 w0 = w · p E · Ly /Lx p E · Ly /Lx 0 t =t· Lx 1 ν0 = ν · p E · Ly Lx

(3)

11 of 14

Francis Chen

User Manual

From the above identities, we then get the rescaling for the derivatives: Lx ∂v 0 ∂v =p 0 ∂x E · Ly /Lx ∂x 0 ∂v Lx ∂v =p ∂y 0 E · Ly /Lx ∂y Lx ∂w0 ∂w =p ∂x0 E · Ly /Lx ∂x 0 Lx ∂w ∂w =p 0 ∂x E · Ly /Lx ∂y

(4)

L2x ∂ 2 v0 ∂2v =p 02 ∂x E · Ly /Lx ∂x2 ∂ 2 v0 L2x ∂2v p = ∂y 02 E · Ly /Lx ∂y 2 ∂ 2 w0 L2x ∂2w =p 02 ∂x E · Ly /Lx ∂x2 ∂ 2 w0 L2x ∂2w p = ∂y 02 E · Ly /Lx ∂y 2

From the above, we can see that the rescaling for velocities and derivatives are totally different, and that is why I use different flags to distinguish different type of outputs. Remember here the initial energy is: initE =

Ly Lx X X

v(i, j)2 + w(i, j)2




(5)

i=0 j=0

There is no

5.3

1 2

here.

Output Format

The output name will be like this: Nx=numOfXGrid Ny=numOfYGrid t=time step ReT=rescaled time step FLAG.txt

Example 6: An example of an output:

Figure 5: Example 6 Here are some explanations for the output data format: 12 of 14

Francis Chen

User Manual

• all the output data are after rescaling. • the first two columns give the rescaled coordinate (x, y). • the last two columns give the values at the coordinate (x, y). The values vary according to the FLAG. See Output Flag. One can refer to the example Output in github repository.

13 of 14

Francis Chen

6

User Manual

Possible Problems • Arrays out of range: the two dimensional arrays in real space is Nx × Ny . In using FFTW, the two dimensional arrays will be converted to 1D, which means the array will be like a[Nx × Ny ]. This may cause the array to be out of range. • Occupying memory: if the grid size is too large, one single output .txt file will be several MB or more. Thus, too many output files will fill the hard disk as a result. To avoid this problem, one can change GENERATE OUTPUT in const.h to a larger one. • Input failure: this should be the problem of wrong input format or wrong input path. For the second problem, one should change the INPUT PATH in const.h to one’s local directory. • Multithreads failure: this should be the constant THREADS in const.h exceeding the maximum threads in one’s machine.

7

Source Code & Documents

All the source code including the FFTW package files can be downloaded from github repository. All the documents can be downloaded from: www.stanford.edu/∼kfchen. The homepage for FFTW is: FFTW Homepage. Any questions regarding the code, please send emails to [email protected].

14 of 14