Commit 766ce3e3 authored by Andreas Nüßing's avatar Andreas Nüßing

[IO] add option to write out the system matrix and the rhs

In order to test some solvers outside of the duneuro-framework, it can be
beneficial to just generate the matrices and the rhss with duneuro, store
them in a file (matrix market) and read the again somewhere else. This
commit adds this option.
Note that this is more of a debug/experimental feature
parent 7c1d6c45
......@@ -12,6 +12,7 @@
#include <dune/common/timer.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixmarket.hh>
#include <dune/pdelab/backend/interface.hh>
#include <dune/pdelab/backend/solver.hh>
......@@ -196,6 +197,7 @@ namespace duneuro
void apply(LS& ls, DV& x, const RV& rightHandSide, const Dune::ParameterTree& config,
DataTree dataTree = DataTree())
{
bool assembledJacobian = false;
Dune::Timer timer(false);
{
std::lock_guard<std::mutex> lock(_jacobian_mutex);
......@@ -209,6 +211,7 @@ namespace duneuro
timer.start();
(*_jacobian) = typename M::field_type(0.0);
_go.jacobian(x, *_jacobian);
assembledJacobian = true;
if (_fixFirstDOF) {
TSSLPDetail::fixFirstDOF(Dune::PDELab::Backend::native(*_jacobian), _fixedDOFEntry);
}
......@@ -238,8 +241,17 @@ namespace duneuro
Dune::Timer solutionTimer;
timer.start();
DV z(_go.trialGridFunctionSpace(), 0.0);
ls.apply(*_jacobian, z, r, config.get<typename RV::ElementType>(
"reduction")); // solver makes right hand side consistent
if (config.get<bool>("only_write_output")) {
if (assembledJacobian) {
std::ofstream matrixStream(config.get<std::string>("matrix_filename"));
Dune::writeMatrixMarket(Dune::PDELab::Backend::native(*_jacobian), matrixStream);
}
std::ofstream rhsStream(config.get<std::string>("rhs_filename"));
Dune::writeMatrixMarket(Dune::PDELab::Backend::native(r), rhsStream);
} else {
ls.apply(*_jacobian, z, r, config.get<typename RV::ElementType>(
"reduction")); // solver makes right hand side consistent
}
timer.stop();
dataTree.set("iterations", ls.result().iterations);
dataTree.set("reduction", ls.result().reduction);
......
......@@ -46,10 +46,12 @@ namespace duneuro
projectedElectrodes.size(), solver_->functionSpace().getGFS().ordering().size());
auto solver_config = config.sub("solver");
typename Traits::DomainDOFVector solution(solver_->functionSpace().getGFS(), 0.0);
auto myConfig = solver_config;
for (std::size_t index = 1; index < projectedElectrodes.size(); ++index) {
myConfig["rhs_filename"] = solver_config["rhs_filename"] + "_" + std::to_string(index);
solve(solverBackend.get(), projectedElectrodes.getProjection(0),
projectedElectrodes.getProjection(index), solution, rightHandSideVector_,
solver_config, dataTree.sub("solver.electrode_" + std::to_string(index)));
myConfig, dataTree.sub("solver.electrode_" + std::to_string(index)));
set_matrix_row(*transferMatrix, index, Dune::PDELab::Backend::native(solution));
}
return transferMatrix;
......@@ -75,10 +77,12 @@ namespace duneuro
tbb::parallel_for(tbb::blocked_range<std::size_t>(1, projectedElectrodes.size(), grainSize),
[&](const tbb::blocked_range<std::size_t>& range) {
auto& mySolution = solution.local();
auto myConfig = solver_config;
for (std::size_t index = range.begin(); index != range.end(); ++index) {
myConfig["rhs_filename"] = solver_config["rhs_filename"] + "_" + std::to_string(index);
solve(solverBackend.local().get(), projectedElectrodes.getProjection(0),
projectedElectrodes.getProjection(index), mySolution,
rightHandSideVector_.local(), solver_config,
rightHandSideVector_.local(), myConfig,
dataTree.sub("solver.electrode_" + std::to_string(index)));
set_matrix_row(*transferMatrix, index,
Dune::PDELab::Backend::native(mySolution));
......
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