Skip to content

Commit

Permalink
Merge pull request #56 from CancerModeling/wagnandr/microcirculation
Browse files Browse the repository at this point in the history 
Wagnandr/microcirculation
  • Loading branch information
@wagnandr
wagnandr authored Mar 24, 2022
2 parents 15ea73f + dc2a523 commit c09df1b
Show file tree
Hide file tree
Showing 38 changed files with 216,929 additions and 96 deletions.
12 changes: 11 additions & 1 deletion apps/breast_geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,13 @@ int main(int argc, char *argv[]) {
auto graph = std::make_shared<mc::GraphStorage>();

mc::EmbeddedGraphReader graph_reader;
graph_reader.append("data/meshes/coarse-network-geometry.json", *graph);
graph_reader.append("data/coarse-breast2-geometry.json", *graph);

// auto inflow_vertices = graph->find_embedded_vertices({ 9.093333333333334, 9.173333333333334, 8.053333333333335 });
std::vector<mc::Point> inflow_points = {
{9.093333333333334, 9.173333333333334, 8.053333333333335},
{2.9333333333333336, 9.973333333333334, 10.933333333333334}};
/*
for (auto &p : inflow_points) {
auto inflow_vertices = graph->find_embedded_vertices(p);
if (inflow_vertices.size() != 1) {
Expand All @@ -45,6 +46,7 @@ int main(int argc, char *argv[]) {
}
inflow_vertices[0]->set_to_inflow_with_fixed_flow(mc::heart_beat_inflow(4.85 / 8.));
}
*/

graph->finalize_bcs();

Expand Down Expand Up @@ -76,11 +78,17 @@ int main(int argc, char *argv[]) {
std::vector<double> p_total_vertex_values;
std::vector<double> p_static_vertex_values;
std::vector<double> c_vertex_values;
std::vector<double> vessel_A0;
std::vector<double> vessel_ids;
std::vector<double> vertex_ids;

mc::fill_with_vessel_A0(MPI_COMM_WORLD, *graph, points, vessel_A0);

// vessels ids do not change, thus we can precalculate them
mc::fill_with_vessel_id(MPI_COMM_WORLD, *graph, points, vessel_ids);

mc::fill_with_vertex_id(MPI_COMM_WORLD, *graph, points, vertex_ids);

// mc::GraphCSVWriter csv_writer(MPI_COMM_WORLD, "output", "data", graph, dof_map, {"Q", "A"});
mc::GraphPVDWriter pvd_writer(MPI_COMM_WORLD, "output", "breast_geometry_solution");

Expand Down Expand Up @@ -110,6 +118,8 @@ int main(int argc, char *argv[]) {
pvd_writer.add_vertex_data("p_total", p_total_vertex_values);
pvd_writer.add_vertex_data("c", c_vertex_values);
pvd_writer.add_vertex_data("vessel_id", vessel_ids);
pvd_writer.add_vertex_data("vertex_ids", vertex_ids);
pvd_writer.add_vertex_data("A0", vessel_A0);
pvd_writer.write(t);
}

Expand Down
74 changes: 37 additions & 37 deletions apps/linear_flow_breast_geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include <petsc.h>
#include <utility>

#include "macrocirculation/0d_boundary_conditions.hpp"
#include "macrocirculation/communication/mpi.hpp"
#include "macrocirculation/dof_map.hpp"
#include "macrocirculation/embedded_graph_reader.hpp"
Expand All @@ -26,39 +27,38 @@ namespace mc = macrocirculation;


int main(int argc, char *argv[]) {
const std::size_t degree = 2;

cxxopts::Options options(argv[0], "Linearized solver for breast geometry");
options.add_options() //
("output-directory", "directory for the output", cxxopts::value<std::string>()->default_value("./output/")) //
("mesh", "filepath to the given mesh", cxxopts::value<std::string>()->default_value("./data/1d-meshes/coarse-breast-geometry.json")) //
("tau", "time step size", cxxopts::value<double>()->default_value("1e-4")) //
("tau-out", "time step size for the output", cxxopts::value<double>()->default_value("1e-2")) //
("t-end", "Endtime for simulation", cxxopts::value<double>()->default_value("10")) //
("h,help", "print usage");
options.allow_unrecognised_options(); // for petsc
auto args = options.parse(argc, argv);
if (args.count("help")) {
std::cout << options.help() << std::endl;
exit(0);
}
if (!args.unmatched().empty()) {
std::cout << "The following arguments were unmatched: " << std::endl;
for (auto &it : args.unmatched())
std::cout << " " << it;
std::cout << "\nAre they part of petsc or a different auxillary library?" << std::endl;
}

const double t_end = args["t-end"].as<double>();
const std::size_t max_iter = 160000000;

const auto tau = args["tau"].as<double>();
const auto tau_out = args["tau-out"].as<double>();

// initialize petsc
CHKERRQ(PetscInitialize(&argc, &argv, nullptr, "solves linear flow problem"));

Eigen::setNbThreads(1);
{
const std::size_t degree = 2;

cxxopts::Options options(argv[0], "Linearized solver for breast geometry");
options.add_options() //
("output-directory", "directory for the output", cxxopts::value<std::string>()->default_value("./output/")) //
("mesh", "filepath to the given mesh", cxxopts::value<std::string>()->default_value("./data/1d-meshes/coarse-breast1-geometry.json")) //
("tau", "time step size", cxxopts::value<double>()->default_value("1e-4")) //
("tau-out", "time step size for the output", cxxopts::value<double>()->default_value("1e-2")) //
("t-end", "Endtime for simulation", cxxopts::value<double>()->default_value("10")) //
("h,help", "print usage");
options.allow_unrecognised_options(); // for petsc
auto args = options.parse(argc, argv);
if (args.count("help")) {
std::cout << options.help() << std::endl;
exit(0);
}
if (!args.unmatched().empty()) {
std::cout << "The following arguments were unmatched: " << std::endl;
for (auto &it : args.unmatched())
std::cout << " " << it;
std::cout << "\nAre they part of petsc or a different auxillary library?" << std::endl;
}

const double t_end = args["t-end"].as<double>();

const auto tau = args["tau"].as<double>();
const auto tau_out = args["tau-out"].as<double>();

std::cout << "rank = " << mc::mpi::rank(PETSC_COMM_WORLD) << std::endl;

const auto output_interval = static_cast<std::size_t>(tau_out / tau);
Expand All @@ -75,8 +75,8 @@ int main(int argc, char *argv[]) {
graph->find_vertex_by_name("bg_135")->set_to_inflow_with_fixed_flow(mc::heart_beat_inflow(0.035));
graph->find_vertex_by_name("bg_141")->set_to_inflow_with_fixed_flow(mc::heart_beat_inflow(0.035));
graph->find_vertex_by_name("bg_119")->set_to_inflow_with_fixed_flow(mc::heart_beat_inflow(0.035));
// mc::set_0d_tree_boundary_conditions(graph, "bg_");
graph_reader.set_boundary_data("./data/meshes/boundary-combined-geometry-linear-part.json", *graph);
mc::set_0d_tree_boundary_conditions(graph, "bg_");
// graph_reader.set_boundary_data("./data/meshes/boundary-combined-geometry-linear-part.json", *graph);

graph->finalize_bcs();

Expand All @@ -86,15 +86,15 @@ int main(int argc, char *argv[]) {
auto dof_map = std::make_shared<mc::DofMap>(graph->num_vertices(), graph->num_edges());
dof_map->create(PETSC_COMM_WORLD, *graph, 2, degree, true);

if (mc::mpi::rank(MPI_COMM_WORLD) == 0)
if (mc::mpi::rank(PETSC_COMM_WORLD) == 0)
std::cout << "num_1d_dof = " << dof_map->num_dof() << std::endl;
std::cout << mc::mpi::rank(MPI_COMM_WORLD) << " owns " << dof_map->num_owned_dofs() << " dof" << std::endl;
std::cout << mc::mpi::rank(PETSC_COMM_WORLD) << " owns " << dof_map->num_owned_dofs() << " dof" << std::endl;

mc::ImplicitLinearFlowSolver solver(PETSC_COMM_WORLD, graph, dof_map, degree);
solver.setup(tau);
solver.use_named_solver("ilf_");

mc::GraphPVDWriter writer(MPI_COMM_WORLD, args["output-directory"].as<std::string>(), "breast_geometry_linearized");
mc::GraphPVDWriter writer(PETSC_COMM_WORLD, args["output-directory"].as<std::string>(), "breast_geometry_linearized");

double t = 0;
const auto t_max_idx = static_cast<size_t>(std::ceil(t_end / tau));
Expand All @@ -119,7 +119,7 @@ int main(int argc, char *argv[]) {
writer.write(t);

int num = 0;
for (auto &v_id : graph->get_active_vertex_ids(mc::mpi::rank(MPI_COMM_WORLD))) {
for (auto &v_id : graph->get_active_vertex_ids(mc::mpi::rank(PETSC_COMM_WORLD))) {
auto &vertex = *graph->get_vertex(v_id);
if (vertex.is_vessel_tree_outflow()) {
const auto &vertex_dof_map = dof_map->get_local_dof_map(vertex);
Expand All @@ -128,7 +128,7 @@ int main(int argc, char *argv[]) {
std::cout << num << std::endl;
num += 1;

std::cout << "rank = " << mc::mpi::rank(MPI_COMM_WORLD) << std::endl;
std::cout << "rank = " << mc::mpi::rank(PETSC_COMM_WORLD) << std::endl;
std::cout << "vertex = " << vertex.get_name() << "\n";
std::cout << " p = ";
for (size_t k = 0; k < vertex_dofs.size(); k += 1)
Expand Down
14 changes: 13 additions & 1 deletion apps/nonlinear_1d_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,6 @@ int main(int argc, char *argv[]) {
}
}

const auto begin_t = std::chrono::steady_clock::now();
double t = 0;
double t_transport = 0.0;

Expand Down Expand Up @@ -187,8 +186,19 @@ int main(int argc, char *argv[]) {

write_output();

double flow_solution_time = 0;
size_t num_iteration = 0;

const auto begin_t = std::chrono::steady_clock::now();

for (std::size_t it = 0; it < max_iter; it += 1) {
auto start = std::chrono::high_resolution_clock::now();
flow_solver->solve(tau, t);
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = (end - start);
flow_solution_time += elapsed.count();
num_iteration += 1;

t += tau;
if (it % update_interval_transport == 0) {
// std::cout << "t_transport = " << t_transport << std::endl;
Expand All @@ -214,6 +224,8 @@ int main(int argc, char *argv[]) {
const auto end_t = std::chrono::steady_clock::now();
const auto elapsed_ms = std::chrono::duration_cast<std::chrono::microseconds>(end_t - begin_t).count();
std::cout << "time = " << elapsed_ms * 1e-6 << " s" << std::endl;

std::cout << "total time flow solver = " << flow_solution_time << ", average = " << flow_solution_time / num_iteration << ", iteration = " << num_iteration << std::endl;
}

CHKERRQ(PetscFinalize());
Expand Down
Loading

0 comments on commit c09df1b

Please sign in to comment.