Skip to content

core: implement scheme as probably wrong forward-euler #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 9 commits into
base: no-bodies
Choose a base branch
from
Open

core: implement scheme as probably wrong forward-euler #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
aa468e7
core: implement scheme as probably wrong forward-euler
blackwer Jan 31, 2023
2fbf8fe
core: add configurational forces to fiber->shell interaction
blackwer Feb 1, 2023
57d0a0d
core: add tension as optional fiber input
blackwer Feb 2, 2023
da33cf2
core: current best guess for what the intended RHS is for fibers...
blackwer Feb 2, 2023
50c0496
core: remove mixing due to reapplication of BC
blackwer Feb 2, 2023
e4d54fb
core: get single fiber + periphery working!
blackwer Feb 2, 2023
7e3f497
core: fix lag in calculating positions that broke self interaction
blackwer Feb 2, 2023
49d54e3
core: change sign on velocities a bit
blackwer Feb 3, 2023
d59f132
core: add motor forces back in
blackwer Feb 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 6 additions & 3 deletions include/fiber.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,9 @@ class Fiber {
x_ = Eigen::MatrixXd::Zero(3, n_nodes_);
x_.row(0) = Eigen::ArrayXd::LinSpaced(n_nodes_, 0, 1.0).transpose();
}
if (tension_.size() != n_nodes_) {
tension_ = Eigen::VectorXd::Zero(n_nodes_);
}

xs_.resize(3, n_nodes_);
xss_.resize(3, n_nodes_);
Expand All @@ -132,7 +135,7 @@ class Fiber {
}


Eigen::VectorXd matvec(VectorRef x, MatrixRef v) const;
Eigen::VectorXd matvec(VectorRef x) const;
void update_preconditioner();
void update_force_operator();
void update_boundary_conditions(Periphery &shell);
Expand Down Expand Up @@ -213,8 +216,8 @@ class FiberContainer {
Eigen::MatrixXd generate_constant_force() const;
const Eigen::MatrixXd &get_local_node_positions() const { return r_fib_local_; };
Eigen::VectorXd get_RHS() const;
Eigen::MatrixXd flow(const MatrixRef &r_trg, const MatrixRef &forces, double eta, bool subtract_self = true) const;
Eigen::VectorXd matvec(VectorRef &x_all, MatrixRef &v_fib) const;
Eigen::MatrixXd flow(const MatrixRef &r_trg, const MatrixRef &forces, double eta, bool subtract_self) const;
Eigen::VectorXd matvec(VectorRef &x_all) const;
Eigen::MatrixXd apply_fiber_force(VectorRef &x_all) const;
Eigen::VectorXd apply_preconditioner(VectorRef &x_all) const;

Expand Down
2 changes: 1 addition & 1 deletion include/periphery.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ class Periphery {
Eigen::VectorXd get_RHS() const { return RHS_; };

Eigen::VectorXd apply_preconditioner(VectorRef &x) const;
Eigen::VectorXd matvec(VectorRef &x_local, MatrixRef &v_local) const;
Eigen::VectorXd matvec(VectorRef &x_local) const;

/// pointer to FMM object (pointer to avoid constructing object with empty Periphery)
kernels::Evaluator stresslet_kernel_;
Expand Down
6 changes: 3 additions & 3 deletions include/system.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ Periphery *get_shell();
PointSourceContainer *get_point_source_container();
toml::value *get_param_table();

std::tuple<int, int> get_local_solution_sizes();
std::tuple<int> get_local_solution_sizes();
Eigen::VectorXd apply_preconditioner(VectorRef &x);
Eigen::VectorXd apply_matvec(VectorRef &x);
void prep_state_for_solver();
Expand All @@ -43,8 +43,8 @@ Eigen::VectorXd get_fiber_RHS();
Eigen::VectorXd get_shell_RHS();
struct properties_t &get_properties();
Eigen::VectorXd &get_curr_solution();
std::tuple<VectorMap, VectorMap> get_solution_maps(double *x);
std::tuple<CVectorMap, CVectorMap> get_solution_maps(const double *x);
std::tuple<VectorMap> get_solution_maps(double *x);
std::tuple<CVectorMap> get_solution_maps(const double *x);
Eigen::MatrixXd velocity_at_targets(MatrixRef &r_trg);

}; // namespace System
Expand Down
34 changes: 4 additions & 30 deletions src/core/fiber.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,10 @@ const std::string Fiber::BC_name[] = {"Force", "Torque", "Velocity", "AngularVel
Fiber::Fiber(toml::value &fiber_table, double eta) {
std::vector<double> x_array = toml::find<std::vector<double>>(fiber_table, "x");
n_nodes_ = x_array.size() / 3;
std::vector<double> tension_array = toml::find_or(fiber_table, "tension", std::vector<double>{});
x_ = Eigen::Map<Eigen::MatrixXd>(x_array.data(), 3, n_nodes_);
if (tension_array.size())
tension_ = Eigen::Map<Eigen::VectorXd>(tension_array.data(), n_nodes_);

init();

Expand Down Expand Up @@ -258,36 +261,7 @@ void Fiber::update_RHS(double dt, MatrixRef &flow, MatrixRef &f_external) {
}
}

VectorXd Fiber::matvec(VectorRef x, MatrixRef v) const {
auto &mats = matrices_.at(n_nodes_);
const int np = n_nodes_;
const int bc_start_i = 4 * np - 14;
MatrixXd D_1 = mats.D_1_0 * std::pow(2.0 / length_, 1);
MatrixXd xsDs = (D_1.array().colwise() * xs_.row(0).transpose().array()).transpose();
MatrixXd ysDs = (D_1.array().colwise() * xs_.row(1).transpose().array()).transpose();
MatrixXd zsDs = (D_1.array().colwise() * xs_.row(2).transpose().array()).transpose();

VectorXd vT = VectorXd(np * 4);
VectorXd v_x = v.row(0).transpose();
VectorXd v_y = v.row(1).transpose();
VectorXd v_z = v.row(2).transpose();

vT.segment(0 * np, np) = v_x;
vT.segment(1 * np, np) = v_y;
vT.segment(2 * np, np) = v_z;
vT.segment(3 * np, np) = xsDs * v_x + ysDs * v_y + zsDs * v_z;

VectorXd vT_in = VectorXd::Zero(4 * np);
vT_in.segment(0, bc_start_i) = mats.P_downsample_bc * vT;

VectorXd xs_vT = VectorXd::Zero(4 * np); // from attachments
const int minus_node = 0;

// FIXME: Does this imply always having velocity boundary conditions?
xs_vT(bc_start_i + 3) = v.col(minus_node).dot(xs_.col(minus_node));

return A_ * x - vT_in + xs_vT;
}
VectorXd Fiber::matvec(VectorRef x) const { return A_ * x; }

/// @brief Calculate the force operator cache variable
///
Expand Down
4 changes: 2 additions & 2 deletions src/core/fiber_container.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,14 +71,14 @@ VectorXd FiberContainer::apply_preconditioner(VectorRef &x_all) const {
return y;
}

VectorXd FiberContainer::matvec(VectorRef &x_all, MatrixRef &v_fib) const {
VectorXd FiberContainer::matvec(VectorRef &x_all) const {
VectorXd res = VectorXd::Zero(get_local_solution_size());

size_t offset = 0;
int i_fib = 0;
for (const auto &fib : *this) {
const int np = fib.n_nodes_;
res.segment(offset, 4 * np) = fib.matvec(x_all.segment(offset, 4 * np), v_fib.block(0, i_fib, 3, np));
res.segment(offset, 4 * np) = fib.matvec(x_all.segment(offset, 4 * np));

i_fib++;
offset += 4 * np;
Expand Down
5 changes: 2 additions & 3 deletions src/core/periphery.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,15 +28,14 @@ Eigen::VectorXd Periphery::apply_preconditioner(VectorRef &x_local) const {
/// @param[in] x_local [3 * n_nodes_local] vector of 'x' local to this rank
/// @param[in] v_local [3 * n_nodes_local] vector of velocities 'v' local to this rank
/// @return [3 * n_nodes_local] vector of A * x_local
Eigen::VectorXd Periphery::matvec(VectorRef &x_local, MatrixRef &v_local) const {
Eigen::VectorXd Periphery::matvec(VectorRef &x_local) const {
if (!n_nodes_global_)
return Eigen::VectorXd();
assert(x_local.size() == get_local_solution_size());
assert(v_local.size() == get_local_solution_size());
Eigen::VectorXd x_shell(3 * n_nodes_global_);
MPI_Allgatherv(x_local.data(), node_counts_[world_rank_], MPI_DOUBLE, x_shell.data(), node_counts_.data(),
node_displs_.data(), MPI_DOUBLE, MPI_COMM_WORLD);
return stresslet_plus_complementary_ * x_shell + CVectorMap(v_local.data(), v_local.size());
return stresslet_plus_complementary_ * x_shell;
}

/// @brief Calculate velocity at target coordinates due to the periphery
Expand Down
12 changes: 5 additions & 7 deletions src/core/solver_hydro.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ P_inv_hydro::P_inv_hydro(const Teuchos::RCP<const Teuchos::Comm<int>> comm) : co
TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), std::invalid_argument,
"P_inv_hydro constructor: The input Teuchos::Comm object must be nonnull.");
const global_ordinal_type indexBase = 0;
const auto [fiber_sol_size, shell_sol_size] = System::get_local_solution_sizes();
const int local_size = fiber_sol_size + shell_sol_size;
const auto [fiber_sol_size] = System::get_local_solution_sizes();
const int local_size = fiber_sol_size;

// Construct a map for our block row distribution
opMap_ = rcp(new map_type(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), local_size, indexBase, comm));
Expand All @@ -32,8 +32,8 @@ A_fiber_hydro::A_fiber_hydro(const Teuchos::RCP<const Teuchos::Comm<int>> comm)
TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), std::invalid_argument,
"A_fiber_hydro constructor: The input Comm object must be nonnull.");
const global_ordinal_type indexBase = 0;
const auto [fiber_sol_size, shell_sol_size] = System::get_local_solution_sizes();
const int local_size = fiber_sol_size + shell_sol_size;
const auto [fiber_sol_size] = System::get_local_solution_sizes();
const int local_size = fiber_sol_size;

// Construct a map for our block row distribution
opMap_ = rcp(new map_type(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), local_size, indexBase, comm));
Expand All @@ -49,13 +49,11 @@ void A_fiber_hydro::apply(const MV &X, MV &Y, Teuchos::ETransp mode, scalar_type

template <>
void Solver<P_inv_hydro, A_fiber_hydro>::set_RHS() {
const auto [fib_sol_size, shell_sol_size] = System::get_local_solution_sizes();
const auto [fib_sol_size] = System::get_local_solution_sizes();
VectorMap RHS_fib(RHS_->getDataNonConst(0).getRawPtr(), fib_sol_size);
VectorMap RHS_shell(RHS_->getDataNonConst(0).getRawPtr() + fib_sol_size, shell_sol_size);

// Initialize GMRES RHS vector
RHS_fib = System::get_fiber_RHS();
RHS_shell = System::get_shell_RHS();
}

template <>
Expand Down
Loading