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

fix sign at boundary

parent e5433fbb
Pipeline #195055 failed with stages
in 23 seconds
This diff is collapsed.
......@@ -95,14 +95,13 @@ namespace finite_differences {
m(0, i-1) = - rxx; // lower
m(1, i-1) = rz + 2. * rxx - f(i); // diagonal
m(2, i-1) = - rxx; // upper;
r(i-1) = (rz - 2. * rxx + fp(i)) * up(i)
+ rxx * up(i-1)
+ rxx * up(i+1);
r(i-1) = (rz + fp(i)) * up(i)
+ rxx * ( up(i-1) - 2 * up(i) + up(i+1));
}
// boundary
r(0) += rxx * u(0);
r(n-1) += rxx * u(n+1);
r(0) -= rxx * u(0);
r(n-1) -= rxx * u(n+1);
// initialize writeable copy
array_1D u_new {u};
......@@ -200,15 +199,14 @@ namespace finite_differences {
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);
r(ix-1) = (rz + fp(ix, iy)) * up(ix, iy)
+ ryy * ( up(ix, iy - 1) - 2.0 * up(ix, iy) + up(ix, iy + 1));
}
// boundary conditions
// TODO(Leon): use interpolated values?
r(0) += rxx * u(0, iy);
r(nx-1) += rxx * u(nx+1, iy);
r(0) -= rxx * u(0, iy);
r(nx-1) -= rxx * u(nx+1, iy);
uhalf.col(iy).segment(1, nx) = algebra::tridiagonal(std::move(m), std::move(r));
}
......@@ -229,13 +227,12 @@ namespace finite_differences {
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(iy-1) = (rz + fp(ix, iy)) * uhalf(ix, iy)
+ rxx * (uhalf(ix-1, iy) - 2.0 * uhalf(ix, iy) + uhalf(ix+1, iy));
}
r(0) += ryy * u(ix, 0);
r(ny-1) += ryy * u(ix, ny+1);
r(0) -= ryy * u(ix, 0);
r(ny-1) -= ryy * u(ix, ny+1);
u_new_trans.col(ix).segment(1, ny) = algebra::tridiagonal(std::move(m), std::move(r));
}
......
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