Commit 89e72c21 authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

use 2+1d solver with constant coefficients

parent 1fd06c6a
This diff is collapsed.
......@@ -203,15 +203,15 @@ class Solver3d():
def _compute_rz(self, Az):
return Az * self._ones
return Az * (1+0j)
def _compute_rxx(self, Axx):
return -Axx * self.dz / 2. / self.dx**2 * self._ones
return -Axx * self.dz / 2. / self.dx**2 * (1+0j)
def _compute_ryy(self, Ayy):
return -Ayy * self.dz / 2. / self.dy**2 * self._ones
return -Ayy * self.dz / 2. / self.dy**2 * (1+0j)
def _compute_f(self, F):
......@@ -235,7 +235,7 @@ class Solver3d():
u[:,-1] = boundary[3]
# the C++ functions use Eigen3 internally, which requires C-Major data storage
u_T = _solver.step2d_A0F(self.rz.T, self.rxx.T, self.ryy.T, fp.T, self.f.T, up.T, u.T)
u_T = _solver.step2d_A0Fs(self.rz, self.rxx, self.ryy, fp.T, self.f.T, up.T, u.T)
self.u = u_T.T
return self.u
......
......@@ -239,4 +239,94 @@ namespace finite_differences{
return u_new;
}
// Scalar Version
// 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,
const complex ryy,
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) {
// 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?
// 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;
r(ix-1) = (rz - 2. * ryy + fp(ix, iy)) * up(ix, iy)
+ ryy * up(ix, iy - 1)
+ ryy * up(ix, iy + 1);
}
// boundary conditions
// TODO: 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);
}
// 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
#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);
}
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);
}
array_2D u_new {u_new_trans.transpose()};
return u_new;
}
}
......@@ -30,6 +30,17 @@ PYBIND11_MODULE(_solver, m){
py::arg("up").noconvert(),
py::arg("u").noconvert()
);
m.def("step2d_A0Fs",
&finite_differences::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()
);
m.def("step1d_AAF",
&finite_differences::step1d_AAF,
......
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