Commit 29123a80 authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

fix non-zero boundary conditions for 3d fd propagator

the boundary conditions for the 3d fd propagator were always 0, even if specified otherwise.
parent 4140124d
......@@ -181,6 +181,15 @@ namespace finite_differences {
array_2D uhalf = array_2D::Zero(u.rows(), u.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));
// 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));
#pragma omp parallel for
for (size_t iy = 1; iy <= ny; ++iy) {
algebra::TriMatrix<scalar> m = algebra::TriMatrix<scalar>::Zero(3, nx);
......@@ -194,9 +203,9 @@ namespace finite_differences {
r(ix-1) = ay * up(ix, iy-1) + (2.0 - 2.0 * ay + 0.5 * fp(ix, iy)) * up(ix, iy) + ay * up(ix, iy+1);
}
// boundary conditions
r(0) += ax * u(0, iy);
r(nx-1) += ax * u(nx+1, iy);
// boundary conditions - interpolate for half-step
r(0) += ax * uhalf(0, iy);
r(nx-1) += ax * uhalf(nx+1, iy);
uhalf.col(iy).segment(1, nx) = algebra::tridiagonal(std::move(m), std::move(r));
}
......@@ -209,6 +218,16 @@ namespace finite_differences {
array_2D u_new_trans = array_2D::Zero(u.cols(), u.rows()); // transposed shape
// 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);
#pragma omp parallel for
for (size_t ix = 1; ix <= nx; ++ix) {
algebra::TriMatrix<scalar> m = algebra::TriMatrix<scalar>::Zero(3, ny);
......
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