Commit ab7a4346 authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

remove unneccesary initialization; change interface for boundary values

parent 6c843169
......@@ -15,7 +15,6 @@ class Solver2d:
self._nx = u0.shape[-1]
self._ones = np.ones((self._nx,), dtype=self.dtype)
self._boundary_slice = np.zeros_like(self._ones)
self.u = u0 * self._ones
......@@ -38,13 +37,7 @@ class Solver2d:
self.f = self._compute_f(F)
u = self._boundary_slice
# set boundary values
u[0] = boundary[0]
u[-1] = boundary[1]
self.u = _solver.step1d_A0F(self.a, fp, self.f, up, u)
self.u = _solver.step1d_A0F(self.a, fp, self.f, up, boundary[0], boundary[1])
return self.u
......@@ -61,8 +54,6 @@ class Solver2dSym:
self._nx = u0.shape[-1]
self._ones = np.ones((self._nx,), dtype=self.dtype)
self._boundary_slice = np.zeros_like(self._ones)
self.u = u0 * self._ones
self.dx = dx
......@@ -84,12 +75,7 @@ class Solver2dSym:
self.f = self._compute_f(F)
u = self._boundary_slice
# set boundary values
u[-1] = boundary[0]
self.u = _solver.step1d_AAF_sym(self.a, fp, self.f, up, u)
self.u = _solver.step1d_AAF_sym(self.a, fp, self.f, up, boundary[0])
return self.u
......@@ -109,7 +95,8 @@ class Solver3d:
self._nx = u0.shape[-1]
self._ny = u0.shape[-2]
self._ones = np.ones(u0.shape, dtype=self.dtype)
self._boundary_slice = np.zeros_like(self._ones)
self._ones_bx = np.ones(u0.shape[-1], dtype=self.dtype)
self._ones_by = np.ones(u0.shape[-2], dtype=self.dtype)
self.u = u0 * self._ones
......@@ -138,16 +125,23 @@ class Solver3d:
self.f = self._compute_f(F)
u = self._boundary_slice
# set boundary values
u[0, :] = boundary[0]
u[-1, :] = boundary[1]
u[:, 0] = boundary[2]
u[:, -1] = boundary[3]
bxl = boundary[0] * self._ones_by
bxr = boundary[1] * self._ones_by
byl = boundary[2] * self._ones_bx
byr = boundary[3] * self._ones_bx
# the C++ functions use Eigen3 internally, which requires Column-Major data storage
u_T = _solver.step2d_A0F(self.ax, self.ay, fp.T, self.f.T, up.T, u.T)
u_T = _solver.step2d_A0F(
self.ax,
self.ay,
fp.T,
self.f.T,
up.T,
bxl,
bxr,
byl,
byr,
)
self.u = u_T.T
return self.u
......@@ -77,18 +77,19 @@ namespace finite_differences {
const Eigen::Ref<const array_1D> fp,
const Eigen::Ref<const array_1D> f,
const Eigen::Ref<const array_1D> up,
const Eigen::Ref<const array_1D> u) {
const scalar bl,
const scalar br ) {
// _p : i
// _ : i+1
size_t n = u.size()-2;
size_t n = up.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);
algebra::TriMatrix<scalar> m (3, n);
array_1D r (n, 1);
// Note that matrix indices go from 1 to u.size()-1, so that indices are off by 1
// Note that matrix indices go from 1 to up.size()-1, so that indices are off by 1
for (size_t i = 1; i <= n; ++i) {
m(0, i-1) = - a; // lower
......@@ -97,12 +98,20 @@ namespace finite_differences {
r(i-1) = a * up(i-1) + (2.0 - 2.0 * a + fp(i)) * up(i) + a * up(i+1);
}
m(0,0) = 0;
m(2,n-1) = 0;
// boundary
r(0) += a * u(0);
r(n-1) += a * u(n+1);
r(0) += a * bl;
r(n-1) += a * br;
// initialize writeable copy
array_1D u_new {u};
array_1D u_new (up.size(), 1);
// insert boundary values
u_new(0) = bl;
u_new(n+1) = br;
u_new.segment(1, n) = algebra::tridiagonal(std::move(m), std::move(r));
return u_new;
......@@ -113,16 +122,16 @@ namespace finite_differences {
const Eigen::Ref<const array_1D> fp,
const Eigen::Ref<const array_1D> f,
const Eigen::Ref<const array_1D> up,
const Eigen::Ref<const array_1D> u) {
const scalar br) {
// _p : i
// _ : i+1
size_t n = u.size()-1;
size_t n = up.size()-1;
// 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);
algebra::TriMatrix<scalar> m (3, n);
array_1D r (n, 1);
for (size_t i = 1; i < n; ++i) {
complex rhix = 1.0 / static_cast<complex>(i);
......@@ -144,10 +153,13 @@ namespace finite_differences {
// right boundary
complex rxnh = 1.0 / static_cast<complex>(n-1);
r(n-1) += (1.0 + 0.5 * rxnh) * a * u(n);
r(n-1) += (1.0 + 0.5 * rxnh) * a * br;
m(2, n-1) = 0;
// initialize writeable copy
array_1D u_new {u};
array_1D u_new (up.size());
u_new(n) = br;
u_new.segment(0, n) = algebra::tridiagonal(std::move(m), std::move(r));
return u_new;
......@@ -155,9 +167,9 @@ namespace finite_differences {
/* 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_A0F(
const complex ax,
......@@ -165,35 +177,37 @@ namespace finite_differences {
const Eigen::Ref<const array_2D> &fp,
const Eigen::Ref<const array_2D> &f,
const Eigen::Ref<const array_2D> &up,
const Eigen::Ref<const array_2D> &u) {
const Eigen::Ref<const array_1D> &bxl,
const Eigen::Ref<const array_1D> &bxr,
const Eigen::Ref<const array_1D> &byl,
const Eigen::Ref<const array_1D> &byr) {
// .p : i
// . : i+1
// 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;
size_t nx = up.rows() - 2;
size_t ny = up.cols() - 2;
// first halfstep
// x at i + 1/2
// y at i
array_2D uhalf = array_2D::Zero(u.rows(), u.cols());
array_2D uhalf (up.rows(), up.cols());
// interpolate boundary values
// x
uhalf.row(0) = 0.5 * (up.row(0) + u.row(0));
uhalf.row(nx+1) = 0.5 * (up.row(nx+1) + u.row(nx+1));
// WARNING: from my understand this should work without the tranpose(). it compiles fine but produces NaNs...
array_1D bxl_half = 0.5 * bxl + 0.5 * up.row(0).transpose();
array_1D bxr_half = 0.5 * bxr + 0.5 * up.row(nx+1).transpose();
array_1D byl_half = 0.5 * byl + 0.5 * up.col(0);
array_1D byr_half = 0.5 * byr + 0.5 * up.col(ny+1);
// y
uhalf.col(0) = 0.5 * (up.col(0) + u.col(0));
uhalf.col(ny+1) = 0.5 * (up.col(ny+1) + u.col(ny+1));
// fill in boundary values. x is done in the loop
uhalf.col(0) = byl_half;
uhalf.col(ny+1) = byr_half;
#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);
algebra::TriMatrix<scalar> m (3, nx);
array_1D r (nx, 1);
for (size_t ix = 1; ix <= nx; ++ix) {
m(0, ix-1) = -ax;
......@@ -204,8 +218,11 @@ namespace finite_differences {
}
// boundary conditions - interpolate for half-step
r(0) += ax * uhalf(0, iy);
r(nx-1) += ax * uhalf(nx+1, iy);
r(0) += ax * bxl_half(iy);
r(nx-1) += ax * bxr_half(iy);
uhalf(0, iy) = bxr_half(iy);
uhalf(nx+1, iy) = bxl_half(iy);
uhalf.col(iy).segment(1, nx) = algebra::tridiagonal(std::move(m), std::move(r));
}
......@@ -214,24 +231,18 @@ namespace finite_differences {
// second halfstep
// note: internal storage is column-major such that column-wise access would be much faster
// TODO(Leon): transpose up, fp, and f
array_2D u_new_trans = array_2D::Zero(u.cols(), u.rows()); // transposed shape
// TODO(Leon): transpose uhalf, fp, and f
// populate boundary values
// x
u_new_trans.col(0) = u.row(0);
u_new_trans.col(nx+1) = u.row(nx+1);
// y
u_new_trans.row(0) = u.col(0);
u_new_trans.row(ny+1) = u.col(ny+1);
array_2D u_new_trans (up.cols(), up.rows()); // transposed shape
// populate boundary values. y is done in the loop
u_new_trans.col(0) = bxl;
u_new_trans.col(nx+1) = bxr;
#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);
algebra::TriMatrix<scalar> m (3, ny);
array_1D r (ny, 1);
for (size_t iy = 1; iy <= ny; ++iy) {
m(0, iy-1) = -ay;
......@@ -242,8 +253,11 @@ namespace finite_differences {
}
// boundary conditions
r(0) += ay * u(ix, 0);
r(ny-1) += ay * u(ix, ny+1);
r(0) += ay * byl(ix);
r(ny-1) += ay * byr(ix);
u_new_trans(0,ix) = byl(ix);
u_new_trans(ny+1,ix) = byr(ix);
u_new_trans.col(ix).segment(1, ny) = algebra::tridiagonal(std::move(m), std::move(r));
}
......
......@@ -13,7 +13,8 @@ PYBIND11_MODULE(_solver, m) {
py::arg("fp").noconvert(),
py::arg("f").noconvert(),
py::arg("up").noconvert(),
py::arg("u").noconvert(),
py::arg("bl").noconvert(),
py::arg("br").noconvert(),
py::return_value_policy::move);
m.def("step2d_A0F",
......@@ -23,7 +24,10 @@ PYBIND11_MODULE(_solver, m) {
py::arg("fp").noconvert(),
py::arg("f").noconvert(),
py::arg("up").noconvert(),
py::arg("u").noconvert(),
py::arg("bxl").noconvert(),
py::arg("bxr").noconvert(),
py::arg("byl").noconvert(),
py::arg("byr").noconvert(),
py::return_value_policy::move);
m.def("step1d_AAF_sym",
......@@ -32,7 +36,7 @@ PYBIND11_MODULE(_solver, m) {
py::arg("fp").noconvert(),
py::arg("f").noconvert(),
py::arg("up").noconvert(),
py::arg("u").noconvert(),
py::arg("bl").noconvert(),
py::return_value_policy::move);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment