diff --git a/src/Compadre_ProblemT.cpp b/src/Compadre_ProblemT.cpp index 7f4e5649f..83452e65f 100644 --- a/src/Compadre_ProblemT.cpp +++ b/src/Compadre_ProblemT.cpp @@ -230,10 +230,6 @@ void ProblemT::solve() { // build solver if (_solver.is_null()) buildSolver(); - //Teuchos::RCP out2 = Teuchos::rcp(new Teuchos::FancyOStream(Teuchos::rcp(new std::ofstream("output.txt")))); - //_solver->getSolution()->describe(*out2,Teuchos::EVerbosityLevel::VERB_EXTREME); - //_A->describe(*out2,Teuchos::EVerbosityLevel::VERB_EXTREME); - // run solver _solver->solve(); @@ -266,6 +262,30 @@ void ProblemT::solve() { } } +void ProblemT::residual() { + + TEUCHOS_TEST_FOR_EXCEPT_MSG(_solver.is_null(), "Can only be called after solve()."); + + for (InteractingFields field_interaction : _field_interactions) { + local_index_type field_one = field_interaction.src_fieldnum; + + local_index_type row_block = _field_to_block_row_map[field_one]; + + _particles->getFieldManager()->updateVectorFromFields(_solver->getSolution(row_block), field_one, _problem_dof_data); + } + + // generate residual, stored in solution + _solver->residual(); + + for (InteractingFields field_interaction : _field_interactions) { + local_index_type field_one = field_interaction.src_fieldnum; + + local_index_type row_block = _field_to_block_row_map[field_one]; + + _particles->getFieldManager()->updateFieldsFromVector(_solver->getSolution(row_block), field_one, _problem_dof_data); + } +} + void ProblemT::buildMaps(local_index_type field_one, local_index_type field_two) { if (field_two<0) field_two = field_one; local_index_type row_block = _field_to_block_row_map[field_one]; diff --git a/src/Compadre_ProblemT.hpp b/src/Compadre_ProblemT.hpp index 9a0ff79e3..cdec7f72d 100644 --- a/src/Compadre_ProblemT.hpp +++ b/src/Compadre_ProblemT.hpp @@ -101,6 +101,7 @@ class ProblemT { void initialize(scalar_type initial_simulation_time = 0); void solve(); + void residual(); private: diff --git a/src/Compadre_SolverT.cpp b/src/Compadre_SolverT.cpp index 52ca50afc..018042d2f 100644 --- a/src/Compadre_SolverT.cpp +++ b/src/Compadre_SolverT.cpp @@ -325,8 +325,6 @@ void SolverT::solve() { //} //printf("SUM IS: %.16f\n", sum); - //auto out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); - ////_A_tpetra[0][0]->describe(*out, Teuchos::VERB_EXTREME); // Constructor from Factory Teuchos::RCP > solver; @@ -343,7 +341,12 @@ void SolverT::solve() { solver->solve(); auto out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); //_A_tpetra[0][0]->describe(*out, Teuchos::VERB_EXTREME); - solver->printTiming(*out, Teuchos::VERB_LOW); + //_A_tpetra[1][0]->describe(*out, Teuchos::VERB_EXTREME); + //_A_single_block->describe(*out, Teuchos::VERB_EXTREME); + //solver->printTiming(*out, Teuchos::VERB_LOW); + //_x_single_block->describe(*out, Teuchos::VERB_EXTREME); + //_b_single_block->describe(*out, Teuchos::VERB_EXTREME); + //_b_tpetra[0]->describe(*out, Teuchos::VERB_EXTREME); } else if (_parameters->get("type")=="iterative") { @@ -440,6 +443,73 @@ void SolverT::solve() { // } } +void SolverT::residual() { + + // can only be called after solve() has been called, because solve() sets up the linear operators + // solution on particles is imported back onto solution vector prior to this call, allowing + // the exact solution to be put on the particles + + if (_consolidate_blocks) { + // move x in blocks into single block, if needed + for (size_t i=0; i<_A_tpetra[0].size(); ++i) { + _x_single_block->doExport(*(_x_tpetra[i]), *(_domain_exporters[i][i]), Tpetra::INSERT); + } + } + + if (_parameters->get("type")=="direct") { + if (_consolidate_blocks) { + _A_single_block->apply(*_x_single_block, *_b_single_block, Teuchos::NO_TRANS, 1.0, -1.0); + // export y->x + _x_single_block->assign(*_b_single_block); + } else { + _A_tpetra[0][0]->apply(*_x_tpetra[0], *_b_tpetra[0], Teuchos::NO_TRANS, 1.0, -1.0); + // export y->x + _x_tpetra[0]->assign(*_b_tpetra[0]); + } + } else { + std::vector x_teko; + if (_consolidate_blocks) { + x_teko.resize(1); + Teuchos::RCP > domain + = Thyra::tpetraVectorSpace(_A_single_block->getDomainMap()); + x_teko[0] = Thyra::tpetraVector(domain, Teuchos::rcp_static_cast(_x_single_block)); + TEUCHOS_TEST_FOR_EXCEPT_MSG(x_teko[0].is_null(), "Cast failed."); + } else { + x_teko.resize(_b_thyra.size()); + for (size_t i=0; i<_b_thyra.size(); i++) { + Teuchos::RCP > domain + = Thyra::tpetraVectorSpace(_A_tpetra[i][i]->getDomainMap()); + x_teko[i] = Thyra::tpetraVector(domain, Teuchos::rcp_static_cast(_x_tpetra[i])); + TEUCHOS_TEST_FOR_EXCEPT_MSG(x_teko[i].is_null(), "Cast failed."); + } + } + + // wrap with Thyra + Teko::MultiVector Thyra_b = Teko::buildBlockedMultiVector(_b_thyra); + Teko::MultiVector Thyra_x = Teko::buildBlockedMultiVector(x_teko); + + _A_thyra->apply(Thyra::NOTRANS, *Thyra_x, Thyra_b.ptr(), 1.0, -1.0); + + // export y->x + if (_consolidate_blocks) { + _x_single_block->assign(*_b_single_block); + } else { + for (size_t i=0; i<_b_thyra.size(); i++) { + _x_tpetra[i]->assign(*_b_tpetra[i]); + } + } + } + + if (_consolidate_blocks) { + // move solution from _x_single_block to _x_tpetra + for (size_t i=0; i<_x_tpetra.size(); i++) { + _x_tpetra[i]->doImport(*_x_single_block, *(_domain_exporters[i][i]), Tpetra::REPLACE); + } + } + + if (_comm->getRank() == 0) std::cout << "WARNING: Over-writing solution variables with the residual." << std::endl; +} + Teuchos::RCP SolverT::getSolution(local_index_type idx) const { Teuchos::RCP return_ptr; @@ -449,6 +519,21 @@ Teuchos::RCP SolverT::getSolution(local_index_type idx) const { } +void SolverT::setSolution(local_index_type idx, Teuchos::RCP vec) { + + auto _x_tpetra_view = _x_tpetra[idx]->getLocalView(); + auto vec_view = vec->getLocalView(); + + TEUCHOS_TEST_FOR_EXCEPT_MSG(_x_tpetra_view.extent(0)!=vec_view.extent(0), "Number of rows differs between vectors."); + TEUCHOS_TEST_FOR_EXCEPT_MSG(_x_tpetra_view.extent(1)!=vec_view.extent(1), "Number of cols differs between vectors."); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,vec_view.extent(0)), KOKKOS_LAMBDA(const int i) { + for (int j=0; j getSolution(local_index_type idx = -1) const; + void setSolution(local_index_type idx, Teuchos::RCP vec); + void solve(); + void residual(); + protected: void amalgamateBlockMatrices(std::vector > > A,