Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
irp
Fresnel
Commits
1e961fd3
Commit
1e961fd3
authored
Feb 17, 2021
by
Leon Merten Lohse
Browse files
add cpplint
parent
fa1223b3
Changes
3
Hide whitespace changes
Inline
Side-by-side
.gitlab-ci.yml
View file @
1e961fd3
...
...
@@ -8,3 +8,9 @@ flake8:
script
:
-
pip install flake8
-
flake8 --max-line-length=120 fresnel/*py
cpplint
:
stage
:
lint
script
:
-
pip install cpplint
-
cpplint --linelength=120 --filter=-legal/copyright,-build/include_order,-build/include_subdir,-runtime/references source/*cpp source/*h
source/finite_difference.h
View file @
1e961fd3
...
...
@@ -2,13 +2,12 @@
#include <complex>
#include <Eigen/Dense>
//#include <iostream>
namespace
algebra
{
namespace
algebra
{
template
<
typename
T
>
using
array1d
=
Eigen
::
Array
<
T
,
Eigen
::
Dynamic
,
1
>
;
template
<
typename
T
>
using
TriMatrix
=
Eigen
::
Array
<
T
,
3
,
Eigen
::
Dynamic
>
;
...
...
@@ -19,18 +18,18 @@ namespace algebra{
*
* Modifies the input coefficients b and r!
*/
template
<
typename
T
>
array1d
<
T
>
tridiagonal
(
TriMatrix
<
T
>
&
m
,
array1d
<
T
>
&
r
)
{
int
n
=
r
.
size
();
auto
a
=
m
.
row
(
0
);
auto
b
=
m
.
row
(
1
);
auto
c
=
m
.
row
(
2
);
array1d
<
T
>
x
(
n
,
1
);
for
(
int
i
=
1
;
i
<
n
;
i
++
)
{
array1d
<
T
>
x
(
n
,
1
);
for
(
int
i
=
1
;
i
<
n
;
i
++
)
{
T
w
=
a
(
i
)
/
b
(
i
-
1
);
b
(
i
)
-=
w
*
c
(
i
-
1
);
r
(
i
)
-=
w
*
r
(
i
-
1
);
...
...
@@ -38,15 +37,15 @@ namespace algebra{
x
(
n
-
1
)
=
r
(
n
-
1
)
/
b
(
n
-
1
);
for
(
int
i
=
n
-
2
;
i
>=
0
;
i
--
)
{
for
(
int
i
=
n
-
2
;
i
>=
0
;
i
--
)
{
x
(
i
)
=
(
r
(
i
)
-
c
(
i
)
*
x
(
i
+
1
))
/
b
(
i
);
}
return
x
;
}
}
}
// namespace algebra
namespace
finite_differences
{
namespace
finite_differences
{
using
real
=
double
;
using
complex
=
std
::
complex
<
real
>
;
...
...
@@ -56,6 +55,7 @@ namespace finite_differences{
using
array_1D
=
Eigen
::
Array
<
scalar
,
Eigen
::
Dynamic
,
1
>
;
using
array_2D
=
Eigen
::
Array
<
scalar
,
Eigen
::
Dynamic
,
Eigen
::
Dynamic
>
;
/*
using shape_t = std::pair<size_t, size_t>;
template <typename T>
...
...
@@ -69,6 +69,7 @@ namespace finite_differences{
throw std::invalid_argument("shape does not match");
dst = src;
}
*/
const
array_1D
step1d_A0F
(
const
Eigen
::
Ref
<
const
array_1D
>
rz
,
...
...
@@ -77,39 +78,38 @@ namespace finite_differences{
const
Eigen
::
Ref
<
const
array_1D
>
f
,
const
Eigen
::
Ref
<
const
array_1D
>
up
,
const
Eigen
::
Ref
<
const
array_1D
>
u
)
{
// _p : i
// _ : i+1
size_t
n
=
u
.
size
()
-
2
;
// _p : i
// _ : i+1
size_t
n
=
u
.
size
()
-
2
;
// setup tridiagonal n x n matrix M = tridiag(a,b,c) and rhs r
algebra
::
TriMatrix
<
scalar
>
m
=
algebra
::
TriMatrix
<
scalar
>::
Zero
(
3
,
n
);
array_1D
r
=
array_1D
::
Zero
(
n
,
1
);
array_1D
r
=
array_1D
::
Zero
(
n
,
1
);
// Note that matrix indices go from 1 to u.size()-1, so that indices are off by 1
for
(
size_t
i
=
1
;
i
<=
n
;
++
i
)
{
m
(
0
,
i
-
1
)
=
-
rxx
(
i
);
// lower
m
(
1
,
i
-
1
)
=
rz
(
i
)
+
2.
*
rxx
(
i
)
-
f
(
i
);
// diagonal
m
(
2
,
i
-
1
)
=
-
rxx
(
i
);
// upper;
r
(
i
-
1
)
=
(
rz
(
i
)
-
2.
*
rxx
(
i
)
+
fp
(
i
)
)
*
up
(
i
)
m
(
0
,
i
-
1
)
=
-
rxx
(
i
);
// lower
m
(
1
,
i
-
1
)
=
rz
(
i
)
+
2.
*
rxx
(
i
)
-
f
(
i
);
// diagonal
m
(
2
,
i
-
1
)
=
-
rxx
(
i
);
// upper;
r
(
i
-
1
)
=
(
rz
(
i
)
-
2.
*
rxx
(
i
)
+
fp
(
i
))
*
up
(
i
)
+
rxx
(
i
)
*
up
(
i
-
1
)
+
rxx
(
i
)
*
up
(
i
+
1
);
}
// boundary
r
(
0
)
+=
rxx
(
1
)
*
u
(
0
);
r
(
n
-
1
)
+=
rxx
(
n
)
*
u
(
n
+
1
);
// initialize writeable copy
array_1D
u_new
{
u
};
u_new
.
segment
(
1
,
n
)
=
algebra
::
tridiagonal
(
m
,
r
);
u_new
.
segment
(
1
,
n
)
=
algebra
::
tridiagonal
(
m
,
r
);
return
u_new
;
}
const
array_1D
step1d_AAF
(
const
Eigen
::
Ref
<
const
array_1D
>
rz
,
const
Eigen
::
Ref
<
const
array_1D
>
rxx
,
...
...
@@ -118,39 +118,37 @@ namespace finite_differences{
const
Eigen
::
Ref
<
const
array_1D
>
f
,
const
Eigen
::
Ref
<
const
array_1D
>
up
,
const
Eigen
::
Ref
<
const
array_1D
>
u
)
{
// _p : i
// _ : i+1
// _p : i
// _ : i+1
size_t
n
=
u
.
size
()
-
2
;
size_t
n
=
u
.
size
()
-
2
;
// setup tridiagonal n x n matrix M = tridiag(a,b,c) and rhs r
algebra
::
TriMatrix
<
scalar
>
m
=
algebra
::
TriMatrix
<
scalar
>::
Zero
(
3
,
n
);
array_1D
r
=
array_1D
::
Zero
(
n
,
1
);
array_1D
r
=
array_1D
::
Zero
(
n
,
1
);
for
(
size_t
i
=
1
;
i
<=
n
;
++
i
)
{
m
(
0
,
i
-
1
)
=
-
rxx
(
i
)
+
rx
(
i
);
// lower
m
(
1
,
i
-
1
)
=
rz
(
i
)
+
2.
*
rxx
(
i
)
-
f
(
i
);
// diagonal
m
(
2
,
i
-
1
)
=
-
rxx
(
i
)
-
rx
(
i
);
// upper;
r
(
i
-
1
)
=
(
rz
(
i
)
-
2.
*
rxx
(
i
)
+
f
(
i
)
)
*
up
(
i
)
m
(
0
,
i
-
1
)
=
-
rxx
(
i
)
+
rx
(
i
);
// lower
m
(
1
,
i
-
1
)
=
rz
(
i
)
+
2.
*
rxx
(
i
)
-
f
(
i
);
// diagonal
m
(
2
,
i
-
1
)
=
-
rxx
(
i
)
-
rx
(
i
);
// upper;
r
(
i
-
1
)
=
(
rz
(
i
)
-
2.
*
rxx
(
i
)
+
f
(
i
))
*
up
(
i
)
+
(
rxx
(
i
)
-
rx
(
i
))
*
up
(
i
-
1
)
+
(
rxx
(
i
)
+
rx
(
i
))
*
up
(
i
+
1
);
}
// boundary
r
(
0
)
+=
(
rxx
(
1
)
-
rx
(
1
))
*
u
(
0
);
r
(
n
-
1
)
+=
(
rxx
(
n
)
+
rx
(
n
))
*
u
(
n
+
1
);
// initialize writeable copy
array_1D
u_new
{
u
};
u_new
.
segment
(
1
,
n
)
=
algebra
::
tridiagonal
(
m
,
r
);
u_new
.
segment
(
1
,
n
)
=
algebra
::
tridiagonal
(
m
,
r
);
return
u_new
;
}
// Two-step alternating-direction Peaceman-Rachford scheme (J. W. Thomas: Numerical Partial Differential Equations)
const
array_2D
step2d_A0F
(
...
...
@@ -161,88 +159,86 @@ namespace finite_differences{
const
Eigen
::
Ref
<
const
array_2D
>
&
f
,
const
Eigen
::
Ref
<
const
array_2D
>
&
up
,
const
Eigen
::
Ref
<
const
array_2D
>
&
u
)
{
// rxx and ryy differ by a factor of 1/2 from the definitions in Fuhse et al. (2005)
// .p : i
// . : i+1
// TODO: interpolate for half-step?
// .p : i
// . : i+1
// TODO
(Leon)
: interpolate for half-step?
// coordinate system:
// numpy (z,y,x), eigen (rows, cols): x <-> cols, y <-> rows
size_t
nx
=
u
.
rows
()
-
2
;
size_t
ny
=
u
.
cols
()
-
2
;
// std::cout << "nx " << u.cols() << ", ny " << u.rows() << std::endl;
// first halfstep
// x at i + 1/2
// y at i
array_2D
uhalf
=
array_2D
::
Zero
(
u
.
rows
(),
u
.
cols
());
#pragma omp parallel for
for
(
size_t
iy
=
1
;
iy
<=
ny
;
++
iy
)
{
algebra
::
TriMatrix
<
scalar
>
m
=
algebra
::
TriMatrix
<
scalar
>::
Zero
(
3
,
nx
);
array_1D
r
=
array_1D
::
Zero
(
nx
,
1
);
for
(
size_t
ix
=
1
;
ix
<=
nx
;
++
ix
)
{
m
(
0
,
ix
-
1
)
=
-
rxx
(
ix
-
1
,
iy
);
m
(
1
,
ix
-
1
)
=
rz
(
ix
,
iy
)
+
2.
*
rxx
(
ix
,
iy
)
-
f
(
ix
,
iy
);
m
(
2
,
ix
-
1
)
=
-
rxx
(
ix
,
iy
);
array_1D
r
=
array_1D
::
Zero
(
nx
,
1
);
for
(
size_t
ix
=
1
;
ix
<=
nx
;
++
ix
)
{
m
(
0
,
ix
-
1
)
=
-
rxx
(
ix
-
1
,
iy
);
m
(
1
,
ix
-
1
)
=
rz
(
ix
,
iy
)
+
2.
*
rxx
(
ix
,
iy
)
-
f
(
ix
,
iy
);
m
(
2
,
ix
-
1
)
=
-
rxx
(
ix
,
iy
);
r
(
ix
-
1
)
=
(
rz
(
ix
,
iy
)
-
2.
*
ryy
(
ix
,
iy
)
+
fp
(
ix
,
iy
))
*
up
(
ix
,
iy
)
+
ryy
(
ix
,
iy
)
*
up
(
ix
,
iy
-
1
)
+
ryy
(
ix
,
iy
)
*
up
(
ix
,
iy
-
1
)
+
ryy
(
ix
,
iy
)
*
up
(
ix
,
iy
+
1
);
}
// boundary conditions
// TODO: use interpolated values?
// TODO
(Leon)
: use interpolated values?
r
(
0
)
+=
rxx
(
1
,
iy
)
*
u
(
0
,
iy
);
r
(
nx
-
1
)
+=
rxx
(
nx
,
iy
)
*
u
(
nx
+
1
,
iy
);
uhalf
.
col
(
iy
).
segment
(
1
,
nx
)
=
algebra
::
tridiagonal
(
m
,
r
);
uhalf
.
col
(
iy
).
segment
(
1
,
nx
)
=
algebra
::
tridiagonal
(
m
,
r
);
}
// second halfstep
// note: internal storage is column-major such that column-wise access would be much faster
array_2D
u_new_trans
=
array_2D
::
Zero
(
u
.
cols
(),
u
.
rows
());
//transposed shape
array_2D
u_new_trans
=
array_2D
::
Zero
(
u
.
cols
(),
u
.
rows
());
//
transposed shape
#pragma omp parallel for
for
(
size_t
ix
=
1
;
ix
<=
nx
;
++
ix
)
{
algebra
::
TriMatrix
<
scalar
>
m
=
algebra
::
TriMatrix
<
scalar
>::
Zero
(
3
,
ny
);
array_1D
r
=
array_1D
::
Zero
(
ny
,
1
);
for
(
size_t
iy
=
1
;
iy
<=
ny
;
++
iy
)
{
m
(
0
,
iy
-
1
)
=
-
ryy
(
ix
,
iy
-
1
);
m
(
1
,
iy
-
1
)
=
rz
(
ix
,
iy
)
+
2.
*
ryy
(
ix
,
iy
)
-
f
(
ix
,
iy
);
m
(
2
,
iy
-
1
)
=
-
ryy
(
ix
,
iy
);
r
(
iy
-
1
)
=
(
rz
(
ix
,
iy
)
-
2.
*
rxx
(
ix
,
iy
)
+
fp
(
ix
,
iy
)
)
*
uhalf
(
ix
,
iy
)
+
rxx
(
ix
,
iy
)
*
uhalf
(
ix
-
1
,
iy
)
+
rxx
(
ix
,
iy
)
*
uhalf
(
ix
+
1
,
iy
);
array_1D
r
=
array_1D
::
Zero
(
ny
,
1
);
for
(
size_t
iy
=
1
;
iy
<=
ny
;
++
iy
)
{
m
(
0
,
iy
-
1
)
=
-
ryy
(
ix
,
iy
-
1
);
m
(
1
,
iy
-
1
)
=
rz
(
ix
,
iy
)
+
2.
*
ryy
(
ix
,
iy
)
-
f
(
ix
,
iy
);
m
(
2
,
iy
-
1
)
=
-
ryy
(
ix
,
iy
);
r
(
iy
-
1
)
=
(
rz
(
ix
,
iy
)
-
2.
*
rxx
(
ix
,
iy
)
+
fp
(
ix
,
iy
))
*
uhalf
(
ix
,
iy
)
+
rxx
(
ix
,
iy
)
*
uhalf
(
ix
-
1
,
iy
)
+
rxx
(
ix
,
iy
)
*
uhalf
(
ix
+
1
,
iy
);
}
r
(
0
)
+=
ryy
(
ix
,
1
)
*
u
(
ix
,
0
);
r
(
ny
-
1
)
+=
ryy
(
ix
,
ny
)
*
u
(
ix
,
ny
+
1
);
u_new_trans
.
col
(
ix
).
segment
(
1
,
ny
)
=
algebra
::
tridiagonal
(
m
,
r
);
u_new_trans
.
col
(
ix
).
segment
(
1
,
ny
)
=
algebra
::
tridiagonal
(
m
,
r
);
}
array_2D
u_new
{
u_new_trans
.
transpose
()};
return
u_new
;
}
// Scalar Version
// Two-step alternating-direction Peaceman-Rachford scheme (J. W. Thomas: Numerical Partial Differential Equations)
/* Two-step alternating-direction Peaceman-Rachford scheme
(J. W. Thomas: Numerical Partial Differential Equations)
*/
const
array_2D
step2d_A0Fs
(
const
complex
rz
,
const
complex
rxx
,
...
...
@@ -251,82 +247,78 @@ namespace finite_differences{
const
Eigen
::
Ref
<
const
array_2D
>
&
f
,
const
Eigen
::
Ref
<
const
array_2D
>
&
up
,
const
Eigen
::
Ref
<
const
array_2D
>
&
u
)
{
// rxx and ryy differ by a factor of 1/2 from the definitions in Fuhse et al. (2005)
// .p : i
// . : i+1
// TODO: interpolate for half-step?
// .p : i
// . : i+1
// TODO
(Leon)
: interpolate for half-step?
// coordinate system:
// numpy (z,y,x), eigen (rows, cols): x <-> cols, y <-> rows
size_t
nx
=
u
.
rows
()
-
2
;
size_t
ny
=
u
.
cols
()
-
2
;
// std::cout << "nx " << u.cols() << ", ny " << u.rows() << std::endl;
// first halfstep
// x at i + 1/2
// y at i
array_2D
uhalf
=
array_2D
::
Zero
(
u
.
rows
(),
u
.
cols
());
#pragma omp parallel for
for
(
size_t
iy
=
1
;
iy
<=
ny
;
++
iy
)
{
algebra
::
TriMatrix
<
scalar
>
m
=
algebra
::
TriMatrix
<
scalar
>::
Zero
(
3
,
nx
);
array_1D
r
=
array_1D
::
Zero
(
nx
,
1
);
for
(
size_t
ix
=
1
;
ix
<=
nx
;
++
ix
)
{
m
(
0
,
ix
-
1
)
=
-
rxx
;
m
(
1
,
ix
-
1
)
=
rz
+
2.
*
rxx
-
f
(
ix
,
iy
);
m
(
2
,
ix
-
1
)
=
-
rxx
;
array_1D
r
=
array_1D
::
Zero
(
nx
,
1
);
for
(
size_t
ix
=
1
;
ix
<=
nx
;
++
ix
)
{
m
(
0
,
ix
-
1
)
=
-
rxx
;
m
(
1
,
ix
-
1
)
=
rz
+
2.
*
rxx
-
f
(
ix
,
iy
);
m
(
2
,
ix
-
1
)
=
-
rxx
;
r
(
ix
-
1
)
=
(
rz
-
2.
*
ryy
+
fp
(
ix
,
iy
))
*
up
(
ix
,
iy
)
+
ryy
*
up
(
ix
,
iy
-
1
)
+
ryy
*
up
(
ix
,
iy
-
1
)
+
ryy
*
up
(
ix
,
iy
+
1
);
}
// boundary conditions
// TODO: use interpolated values?
// TODO
(Leon)
: use interpolated values?
r
(
0
)
+=
rxx
*
u
(
0
,
iy
);
r
(
nx
-
1
)
+=
rxx
*
u
(
nx
+
1
,
iy
);
uhalf
.
col
(
iy
).
segment
(
1
,
nx
)
=
algebra
::
tridiagonal
(
m
,
r
);
uhalf
.
col
(
iy
).
segment
(
1
,
nx
)
=
algebra
::
tridiagonal
(
m
,
r
);
}
// second halfstep
// note: internal storage is column-major such that column-wise access would be much faster
array_2D
u_new_trans
=
array_2D
::
Zero
(
u
.
cols
(),
u
.
rows
());
//transposed shape
array_2D
u_new_trans
=
array_2D
::
Zero
(
u
.
cols
(),
u
.
rows
());
//
transposed shape
#pragma omp parallel for
for
(
size_t
ix
=
1
;
ix
<=
nx
;
++
ix
)
{
algebra
::
TriMatrix
<
scalar
>
m
=
algebra
::
TriMatrix
<
scalar
>::
Zero
(
3
,
ny
);
array_1D
r
=
array_1D
::
Zero
(
ny
,
1
);
for
(
size_t
iy
=
1
;
iy
<=
ny
;
++
iy
)
{
m
(
0
,
iy
-
1
)
=
-
ryy
;
m
(
1
,
iy
-
1
)
=
rz
+
2.
*
ryy
-
f
(
ix
,
iy
);
m
(
2
,
iy
-
1
)
=
-
ryy
;
r
(
iy
-
1
)
=
(
rz
-
2.
*
rxx
+
fp
(
ix
,
iy
)
)
*
uhalf
(
ix
,
iy
)
+
rxx
*
uhalf
(
ix
-
1
,
iy
)
+
rxx
*
uhalf
(
ix
+
1
,
iy
);
array_1D
r
=
array_1D
::
Zero
(
ny
,
1
);
for
(
size_t
iy
=
1
;
iy
<=
ny
;
++
iy
)
{
m
(
0
,
iy
-
1
)
=
-
ryy
;
m
(
1
,
iy
-
1
)
=
rz
+
2.
*
ryy
-
f
(
ix
,
iy
);
m
(
2
,
iy
-
1
)
=
-
ryy
;
r
(
iy
-
1
)
=
(
rz
-
2.
*
rxx
+
fp
(
ix
,
iy
))
*
uhalf
(
ix
,
iy
)
+
rxx
*
uhalf
(
ix
-
1
,
iy
)
+
rxx
*
uhalf
(
ix
+
1
,
iy
);
}
r
(
0
)
+=
ryy
*
u
(
ix
,
0
);
r
(
ny
-
1
)
+=
ryy
*
u
(
ix
,
ny
+
1
);
u_new_trans
.
col
(
ix
).
segment
(
1
,
ny
)
=
algebra
::
tridiagonal
(
m
,
r
);
u_new_trans
.
col
(
ix
).
segment
(
1
,
ny
)
=
algebra
::
tridiagonal
(
m
,
r
);
}
array_2D
u_new
{
u_new_trans
.
transpose
()};
return
u_new
;
}
}
}
// namespace finite_differences
source/solver_py.cpp
View file @
1e961fd3
...
...
@@ -6,50 +6,46 @@
namespace
py
=
pybind11
;
using
namespace
finite_differences
;
namespace
fd
=
finite_differences
;
PYBIND11_MODULE
(
_solver
,
m
){
PYBIND11_MODULE
(
_solver
,
m
)
{
m
.
def
(
"step1d_A0F"
,
&
f
inite_differences
::
step1d_A0F
,
&
f
d
::
step1d_A0F
,
py
::
arg
(
"rz"
).
noconvert
(),
py
::
arg
(
"rxx"
).
noconvert
(),
py
::
arg
(
"fp"
).
noconvert
(),
py
::
arg
(
"f"
).
noconvert
(),
py
::
arg
(
"up"
).
noconvert
(),
py
::
arg
(
"u"
).
noconvert
()
);
py
::
arg
(
"u"
).
noconvert
()
);
m
.
def
(
"step2d_A0F"
,
&
f
inite_differences
::
step2d_A0F
,
&
f
d
::
step2d_A0F
,
py
::
arg
(
"rz"
).
noconvert
(),
py
::
arg
(
"rxx"
).
noconvert
(),
py
::
arg
(
"ryy"
).
noconvert
(),
py
::
arg
(
"fp"
).
noconvert
(),
py
::
arg
(
"f"
).
noconvert
(),
py
::
arg
(
"up"
).
noconvert
(),
py
::
arg
(
"u"
).
noconvert
()
);
py
::
arg
(
"u"
).
noconvert
()
);
m
.
def
(
"step2d_A0Fs"
,
&
f
inite_differences
::
step2d_A0Fs
,
&
f
d
::
step2d_A0Fs
,
py
::
arg
(
"rz"
).
noconvert
(),
py
::
arg
(
"rxx"
).
noconvert
(),
py
::
arg
(
"ryy"
).
noconvert
(),
py
::
arg
(
"fp"
).
noconvert
(),
py
::
arg
(
"f"
).
noconvert
(),
py
::
arg
(
"up"
).
noconvert
(),
py
::
arg
(
"u"
).
noconvert
()
);
py
::
arg
(
"u"
).
noconvert
()
);
m
.
def
(
"step1d_AAF"
,
&
f
inite_differences
::
step1d_AAF
,
&
f
d
::
step1d_AAF
,
py
::
arg
(
"rz"
).
noconvert
(),
py
::
arg
(
"rxx"
).
noconvert
(),
py
::
arg
(
"rx"
).
noconvert
(),
py
::
arg
(
"fp"
).
noconvert
(),
py
::
arg
(
"f"
).
noconvert
(),
py
::
arg
(
"up"
).
noconvert
(),
py
::
arg
(
"u"
).
noconvert
()
);
}
\ No newline at end of file
py
::
arg
(
"u"
).
noconvert
()
);
}
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment