diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp index a5e3d75d2362..60d1c6d87f75 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp @@ -28,7 +28,9 @@ using Teuchos::rcp; #include "Panzer_GatherOrientation.hpp" #include "Panzer_ScatterResidual_BlockedTpetra.hpp" #include "Panzer_GatherSolution_BlockedTpetra.hpp" +#include "Panzer_LOCPair_GlobalEvaluationData.hpp" #include "Panzer_GlobalEvaluationDataContainer.hpp" +#include "Panzer_ParameterList_GlobalEvaluationData.hpp" #include "Panzer_STK_Version.hpp" #include "PanzerAdaptersSTK_config.hpp" @@ -494,6 +496,367 @@ namespace panzer } } + TEUCHOS_UNIT_TEST(block_assembly, scatter_solution_tangent) + { + +#ifdef HAVE_MPI + Teuchos::RCP> tComm = Teuchos::rcp(new Teuchos::MpiComm(MPI_COMM_WORLD)); +#else + NOPE_PANZER_DOESNT_SUPPORT_SERIAL +#endif + + int myRank = tComm->getRank(); + + const std::size_t workset_size = 4; + const std::string fieldName1_q1 = "U"; + const std::string fieldName2_q1 = "V"; + const std::string fieldName_qedge1 = "B"; + const std::size_t numParams = 2; + + Teuchos::RCP mesh = buildMesh(2, 2); + + // build input physics block + Teuchos::RCP basis_q1 = buildBasis(workset_size, "Q1"); + Teuchos::RCP basis_qedge1 = buildBasis(workset_size, "QEdge1"); + + Teuchos::RCP ipb = Teuchos::parameterList(); + testInitialization(ipb); + + const int default_int_order = 1; + std::string eBlockID = "eblock-0_0"; + Teuchos::RCP eqset_factory = Teuchos::rcp(new user_app::MyFactory); + panzer::CellData cellData(workset_size, mesh->getCellTopology("eblock-0_0")); + Teuchos::RCP gd = panzer::createGlobalData(); + Teuchos::RCP physicsBlock = + Teuchos::rcp(new PhysicsBlock(ipb, eBlockID, default_int_order, cellData, eqset_factory, gd, false)); + + Teuchos::RCP> work_sets = panzer_stk::buildWorksets(*mesh, physicsBlock->elementBlockID(), + physicsBlock->getWorksetNeeds()); + TEST_EQUALITY(work_sets->size(), 1); + + // build connection manager and field manager + const Teuchos::RCP conn_manager = Teuchos::rcp(new panzer_stk::STKConnManager(mesh)); + RCP dofManager = Teuchos::rcp(new panzer::BlockedDOFManager(conn_manager, MPI_COMM_WORLD)); + + dofManager->addField(fieldName1_q1, Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis()))); + dofManager->addField(fieldName2_q1, Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis()))); + dofManager->addField(fieldName_qedge1, Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_qedge1->getIntrepid2Basis()))); + + std::vector> fieldOrder(3); + fieldOrder[0].push_back(fieldName1_q1); + fieldOrder[1].push_back(fieldName_qedge1); + fieldOrder[2].push_back(fieldName2_q1); + dofManager->setFieldOrder(fieldOrder); + + // dofManager->setOrientationsRequired(true); + dofManager->buildGlobalUnknowns(); + + // setup linear object factory + ///////////////////////////////////////////////////////////// + Teuchos::RCP bt_lof = Teuchos::rcp(new TpetraBlockedLinObjFactoryType(tComm.getConst(), dofManager)); + Teuchos::RCP> lof = bt_lof; + Teuchos::RCP loc = bt_lof->buildGhostedLinearObjContainer(); + bt_lof->initializeGhostedContainer(LinearObjContainer::X | LinearObjContainer::F, *loc); + loc->initialize(); + + Teuchos::RCP b_loc = Teuchos::rcp_dynamic_cast(loc); + + Teuchos::RCP> p_vec = Teuchos::rcp_dynamic_cast>(b_loc->get_x()); + Thyra::assign(p_vec->getNonconstVectorBlock(0).ptr(), 123.0 + myRank); + Thyra::assign(p_vec->getNonconstVectorBlock(1).ptr(), 456.0 + myRank); + Thyra::assign(p_vec->getNonconstVectorBlock(2).ptr(), 789.0 + myRank); + + std::vector> tangentContainers; + + using LOCPair = panzer::LOCPair_GlobalEvaluationData; + using Teuchos::rcp_dynamic_cast; + + // generate tangent data + for (std::size_t i=0;i(locPair->getGlobalLOC()); + Teuchos::RCP> global_p_vec = Teuchos::rcp_dynamic_cast>(global_bt_loc->get_x()); + Thyra::assign(global_p_vec->getNonconstVectorBlock(0).ptr(), 0.123 + myRank + i); + Thyra::assign(global_p_vec->getNonconstVectorBlock(1).ptr(), 0.456 + myRank + i); + Thyra::assign(global_p_vec->getNonconstVectorBlock(2).ptr(), 0.789 + myRank + i); + + auto ghosted_bt_loc = rcp_dynamic_cast(locPair->getGhostedLOC()); + Teuchos::RCP> ghosted_p_vec = Teuchos::rcp_dynamic_cast>(ghosted_bt_loc->get_x()); + Thyra::assign(ghosted_p_vec->getNonconstVectorBlock(0).ptr(), 0.123 + myRank + i); + Thyra::assign(ghosted_p_vec->getNonconstVectorBlock(1).ptr(), 0.456 + myRank + i); + Thyra::assign(ghosted_p_vec->getNonconstVectorBlock(2).ptr(), 0.789 + myRank + i); + + tangentContainers.push_back(locPair); + } + + // setup field manager, add evaluator under test + ///////////////////////////////////////////////////////////// + + auto fm = Teuchos::rcp(new PHX::FieldManager); + + std::vector derivative_dimensions; + derivative_dimensions.push_back(numParams); + fm->setKokkosExtendedDataTypeDimensions(derivative_dimensions); + + std::string resName = ""; + Teuchos::RCP> names_map = + Teuchos::rcp(new std::map); + names_map->insert(std::make_pair(fieldName1_q1, resName + fieldName1_q1)); + names_map->insert(std::make_pair(fieldName2_q1, resName + fieldName2_q1)); + names_map->insert(std::make_pair(fieldName_qedge1, resName + fieldName_qedge1)); + + // evaluators under test + { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + names->push_back(resName + fieldName1_q1); + names->push_back(resName + fieldName2_q1); + + Teuchos::ParameterList pl; + pl.set("Scatter Name", "ScatterQ1"); + pl.set("Basis", basis_q1.getConst()); + pl.set("Dependent Names", names); + pl.set("Dependent Map", names_map); + + Teuchos::RCP> evaluator = lof->buildScatter(pl); + + TEST_EQUALITY(evaluator->evaluatedFields().size(), 1); + + fm->registerEvaluator(evaluator); + fm->requireField(*evaluator->evaluatedFields()[0]); + } + { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + names->push_back(resName + fieldName_qedge1); + + Teuchos::ParameterList pl; + pl.set("Scatter Name", "ScatterQEdge1"); + pl.set("Basis", basis_qedge1.getConst()); + pl.set("Dependent Names", names); + pl.set("Dependent Map", names_map); + + Teuchos::RCP> evaluator = lof->buildScatter(pl); + + TEST_EQUALITY(evaluator->evaluatedFields().size(), 1); + + fm->registerEvaluator(evaluator); + fm->requireField(*evaluator->evaluatedFields()[0]); + } + + // support evaluators + { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + names->push_back(fieldName1_q1); + names->push_back(fieldName2_q1); + + Teuchos::ParameterList pl; + pl.set("Basis", basis_q1); + pl.set("DOF Names", names); + pl.set("Indexer Names", names); + Teuchos::RCP>> tangent_names = + Teuchos::rcp(new std::vector>(2)); + for (std::size_t i = 0; i < numParams; ++i) + { + std::stringstream ss1, ss2; + ss1 << fieldName1_q1 << " Tangent " << i; + ss2 << fieldName2_q1 << " Tangent " << i; + (*tangent_names)[0].push_back(ss1.str()); + (*tangent_names)[1].push_back(ss2.str()); + } + pl.set("Tangent Names", tangent_names); + + Teuchos::RCP> evaluator = lof->buildGather(pl); + + fm->registerEvaluator(evaluator); + } + for (std::size_t i = 0; i < numParams; ++i) { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + RCP> tangent_names = rcp(new std::vector); + names->push_back(fieldName1_q1); + names->push_back(fieldName2_q1); + { + std::stringstream ss1, ss2; + ss1 << fieldName1_q1 << " Tangent " << i; + ss2 << fieldName2_q1 << " Tangent " << i; + tangent_names->push_back(ss1.str()); + tangent_names->push_back(ss2.str()); + } + + Teuchos::ParameterList pl; + pl.set("Basis", basis_q1); + pl.set("DOF Names", tangent_names); + pl.set("Indexer Names", names); + + std::stringstream ss; + ss << "Tangent Container " << i; + pl.set("Global Data Key", ss.str()); + + Teuchos::RCP> evaluator = + lof->buildGatherTangent(pl); + + fm->registerEvaluator(evaluator); + } + { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + names->push_back(fieldName_qedge1); + + Teuchos::ParameterList pl; + pl.set("Basis", basis_qedge1); + pl.set("DOF Names", names); + pl.set("Indexer Names", names); + Teuchos::RCP>> tangent_names = + Teuchos::rcp(new std::vector>(1)); + for (std::size_t i = 0; i < numParams; ++i) + { + std::stringstream ss; + ss << fieldName_qedge1 << " Tangent " << i; + (*tangent_names)[0].push_back(ss.str()); + } + pl.set("Tangent Names", tangent_names); + + Teuchos::RCP> evaluator = lof->buildGather(pl); + + fm->registerEvaluator(evaluator); + } + for (std::size_t i = 0; i < numParams; ++i) { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + RCP> tangent_names = rcp(new std::vector); + names->push_back(fieldName_qedge1); + { + std::stringstream ss; + ss << fieldName_qedge1 << " Tangent " << i; + tangent_names->push_back(ss.str()); + } + + Teuchos::ParameterList pl; + pl.set("Basis", basis_qedge1); + pl.set("DOF Names", tangent_names); + pl.set("Indexer Names", names); + + std::stringstream ss; + ss << "Tangent Container " << i; + pl.set("Global Data Key", ss.str()); + + Teuchos::RCP> evaluator = + lof->buildGatherTangent(pl); + + fm->registerEvaluator(evaluator); + } + + panzer::Traits::SD sd; + sd.worksets_ = work_sets; + + fm->postRegistrationSetup(sd); + + panzer::Traits::PED ped; + ped.gedc->addDataObject("Solution Gather Container", loc); + ped.gedc->addDataObject("Residual Scatter Container", loc); + for (size_t i=0; iaddDataObject(ss.str(), tangentContainers[i]); + } + std::vector params; + std::vector> paramContainers; + for (std::size_t i = 0; iaddDataObject(ss.str(),paramContainer->getGhostedLOC()); + paramContainers.push_back(paramContainer); + } + Teuchos::RCP activeParams = + Teuchos::rcp(new panzer::ParameterList_GlobalEvaluationData(params)); + ped.gedc->addDataObject("PARAMETER_NAMES",activeParams); + fm->preEvaluate(ped); + + // run tests + ///////////////////////////////////////////////////////////// + + panzer::Workset &workset = (*work_sets)[0]; + workset.alpha = 0.0; + workset.beta = 2.0; // derivatives multiplied by 2 + workset.time = 0.0; + workset.evaluate_transient_terms = false; + + fm->evaluateFields(workset); + fm->postEvaluate(0); + + fm = Teuchos::null; + + // test Tangent fields + Teuchos::ArrayRCP data; + Teuchos::RCP> f_vec = Teuchos::rcp_dynamic_cast>(b_loc->get_f()); + + // check all the residual values. This is kind of crappy test since it simply checks twice the target + // value and the target. Its this way because you add two entries across elements. + + Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(0))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(0)->getLocalNumElements()); + for (size_type i = 0; i < data.size(); i++) + { + double target = 123.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + } + + Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(1))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(1)->getLocalNumElements()); + for (size_type i = 0; i < data.size(); i++) + { + double target = 456.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + } + + Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(2))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(2)->getLocalNumElements()); + for (size_type i = 0; i < data.size(); i++) + { + double target = 789.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + } + + for (std::size_t i=0; i> param_f_vec = + Teuchos::rcp_dynamic_cast>( + Teuchos::rcp_dynamic_cast(paramContainers[i]->getGhostedLOC())->get_f()); + + Teuchos::rcp_dynamic_cast>(param_f_vec->getVectorBlock(0))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(0)->getLocalNumElements()); + for (size_type j = 0; j < data.size(); j++) + { + double target = .123 + myRank + i; + TEST_ASSERT(data[j] == target || data[j] == 2.0 * target); + } + Teuchos::rcp_dynamic_cast>(param_f_vec->getVectorBlock(1))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(1)->getLocalNumElements()); + for (size_type j = 0; j < data.size(); j++) + { + double target = .456 + myRank + i; + TEST_ASSERT(data[j] == target || data[j] == 2.0 * target); + } + Teuchos::rcp_dynamic_cast>(param_f_vec->getVectorBlock(2))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(2)->getLocalNumElements()); + for (size_type j = 0; j < data.size(); j++) + { + double target = .789 + myRank + i; + TEST_ASSERT(data[j] == target || data[j] == 2.0 * target); + } + } + } + Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName) { Teuchos::RCP topo = diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra.hpp index bb0522da313a..382c27d64be0 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra.hpp @@ -51,7 +51,7 @@ class GatherTangent_BlockedTpetra public: GatherTangent_BlockedTpetra(const Teuchos::RCP & indexer) - : gidIndexer_(indexer) {} + : globalIndexer_(indexer) {} GatherTangent_BlockedTpetra(const Teuchos::RCP & indexer, const Teuchos::ParameterList& p); @@ -64,13 +64,13 @@ class GatherTangent_BlockedTpetra void evaluateFields(typename TRAITS::EvalData d); virtual Teuchos::RCP clone(const Teuchos::ParameterList & pl) const - { return Teuchos::rcp(new GatherTangent_BlockedTpetra(gidIndexer_,pl)); } + { return Teuchos::rcp(new GatherTangent_BlockedTpetra(globalIndexer_,pl)); } private: // We always use RealType for gathering as we never compute derivatives for this evaluator - //typedef typename panzer::Traits::RealType ScalarT; - typedef typename EvalT::ScalarT ScalarT; + typedef typename panzer::Traits::RealType ScalarT; + //typedef typename EvalT::ScalarT ScalarT; typedef BlockedTpetraLinearObjContainer ContainerType; typedef Tpetra::Vector VectorType; @@ -82,10 +82,17 @@ class GatherTangent_BlockedTpetra // maps the local (field,element,basis) triplet to a global ID // for scattering - Teuchos::RCP gidIndexer_; + Teuchos::RCP globalIndexer_; std::vector fieldIds_; // field IDs needing mapping + //! Vector of global indexers, one for each field to gather, respectively + std::vector> fieldGlobalIndexers_; + + //! Returns the index to the Thyra ProductVector sub-block. Size + //! of number of fields to gather + std::vector productVectorBlockIndex_; + std::vector< PHX::MDField > gatherFields_; Teuchos::RCP > indexerNames_; @@ -94,6 +101,12 @@ class GatherTangent_BlockedTpetra Teuchos::RCP > blockedContainer_; + //! Local indices for unknowns + PHX::View worksetLIDs_; + + //! Offset into the cell lids for each field + std::vector> fieldOffsets_; + GatherTangent_BlockedTpetra(); }; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra_impl.hpp index adf75295efe8..6673e3bc8553 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherTangent_BlockedTpetra_impl.hpp @@ -18,6 +18,7 @@ #include "Panzer_BlockedDOFManager.hpp" #include "Panzer_PureBasis.hpp" #include "Panzer_TpetraLinearObjFactory.hpp" +#include "Panzer_LOCPair_GlobalEvaluationData.hpp" #include "Panzer_BlockedTpetraLinearObjContainer.hpp" #include "Panzer_GlobalEvaluationDataContainer.hpp" @@ -33,7 +34,7 @@ panzer::GatherTangent_BlockedTpetra:: GatherTangent_BlockedTpetra( const Teuchos::RCP & indexer, const Teuchos::ParameterList& p) - : gidIndexer_(indexer) + : globalIndexer_(indexer) , useTimeDerivativeSolutionVector_(false) , globalDataKey_("Tangent Gather Container") { @@ -50,6 +51,10 @@ GatherTangent_BlockedTpetra( gatherFields_[fd] = PHX::MDField(names[fd],basis->functional); this->addEvaluatedField(gatherFields_[fd]); + // If blockedContainer_ is null, the evalaution is a no-op. In this + // case we need to preserve zero initial value. Do this by not + // sharing. + this->addUnsharedField(gatherFields_[fd].fieldTag().clone()); } if (p.isType("Use Time Derivative Solution Vector")) @@ -64,19 +69,44 @@ GatherTangent_BlockedTpetra( // ********************************************************************** template void panzer::GatherTangent_BlockedTpetra:: -postRegistrationSetup(typename TRAITS::SetupData /* d */, +postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager& /* fm */) { TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_->size()); - fieldIds_.resize(gatherFields_.size()); + const Workset & workset_0 = (*d.worksets_)[0]; + const std::string blockId = this->wda(workset_0).block_id; + fieldIds_.resize(gatherFields_.size()); + fieldOffsets_.resize(gatherFields_.size()); + fieldGlobalIndexers_.resize(gatherFields_.size()); + productVectorBlockIndex_.resize(gatherFields_.size()); + int maxElementBlockGIDCount = -1; for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) { // get field ID from DOF manager const std::string& fieldName = (*indexerNames_)[fd]; - fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName); + const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager + productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum); + fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]]; + fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName); // Field number in the sub-global-indexer + + const std::vector& offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]); + fieldOffsets_[fd] = PHX::View("GatherSolution_BlockedTpetra(Residual):fieldOffsets",offsets.size()); + auto hostFieldOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]); + for(std::size_t i=0; i < offsets.size(); ++i) + hostFieldOffsets(i) = offsets[i]; + Kokkos::deep_copy(fieldOffsets_[fd],hostFieldOffsets); + + maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount); } + // We will use one workset lid view for all fields, but has to be + // sized big enough to hold the largest elementBlockGIDCount in the + // ProductVector. + worksetLIDs_ = PHX::View("GatherSolution_BlockedTpetra(Residual):worksetLIDs", + gatherFields_[0].extent(0), + maxElementBlockGIDCount); + indexerNames_ = Teuchos::null; // Don't need this anymore } @@ -87,7 +117,18 @@ preEvaluate(typename TRAITS::PreEvalData d) { // try to extract linear object container if (d.gedc->containsDataObject(globalDataKey_)) { - blockedContainer_ = Teuchos::rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_),true); + Teuchos::RCP ged = d.gedc->getDataObject(globalDataKey_); + Teuchos::RCP loc_pair = + Teuchos::rcp_dynamic_cast(ged); + + if(loc_pair!=Teuchos::null) { + Teuchos::RCP loc = loc_pair->getGhostedLOC(); + blockedContainer_ = Teuchos::rcp_dynamic_cast(loc,true); + } + + if(blockedContainer_==Teuchos::null) { + blockedContainer_ = Teuchos::rcp_dynamic_cast(ged,true); + } } } @@ -96,73 +137,50 @@ template :: evaluateFields(typename TRAITS::EvalData workset) { - using Teuchos::RCP; - using Teuchos::ArrayRCP; - using Teuchos::ptrFromRef; - using Teuchos::rcp_dynamic_cast; - - using Thyra::VectorBase; - using Thyra::SpmdVectorBase; - using Thyra::ProductVectorBase; - - // If blockedContainer_ was not initialized, then no global evaluation data - // container was set, in which case this evaluator becomes a no-op - if (blockedContainer_ == Teuchos::null) - return; - - Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout)); - out.setShowProcRank(true); - out.setOutputToRootOnly(-1); - - std::vector > GIDs; - std::vector LIDs; - - // for convenience pull out some objects from workset - std::string blockId = this->wda(workset).block_id; - const std::vector & localCellIds = this->wda(workset).cell_local_ids; - - Teuchos::RCP > x; - if (useTimeDerivativeSolutionVector_) - x = rcp_dynamic_cast >(blockedContainer_->get_dxdt()); - else - x = rcp_dynamic_cast >(blockedContainer_->get_x()); - - // gather operation for each cell in workset - for(std::size_t worksetCellIndex=0;worksetCellIndexgetElementGIDsPair(cellLocalId,GIDs,blockId); - - // caculate the local IDs for this element - LIDs.resize(GIDs.size()); - for(std::size_t i=0;i x_map = blockedContainer_->getMapForBlock(GIDs[i].first); - - LIDs[i] = x_map->getLocalElement(GIDs[i].second); + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Thyra::VectorBase; + using Thyra::ProductVectorBase; + + // If blockedContainer_ was not initialized, then no global evaluation data + // container was set, in which case this evaluator becomes a no-op + if (blockedContainer_ == Teuchos::null) return; + + const PHX::View& localCellIds = this->wda(workset).cell_local_ids_k; + + RCP> thyraBlockSolution; + if (useTimeDerivativeSolutionVector_) + thyraBlockSolution = rcp_dynamic_cast>(blockedContainer_->get_dxdt(),true); + else + thyraBlockSolution = rcp_dynamic_cast>(blockedContainer_->get_x(),true); + + // Loop over gathered fields + int currentWorksetLIDSubBlock = -1; + for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) { + // workset LIDs only change for different sub blocks + if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) { + const std::string blockId = this->wda(workset).block_id; + const int num_dofs = fieldGlobalIndexers_[fieldIndex]->getElementBlockGIDCount(blockId); + fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_,num_dofs); + currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex]; + } + + const auto& tpetraSolution = *((rcp_dynamic_cast>(thyraBlockSolution->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),true))->getTpetraVector()); + const auto& kokkosSolution = tpetraSolution.getLocalViewDevice(Tpetra::Access::ReadOnly); + + // Class data fields for lambda capture + const auto& fieldOffsets = fieldOffsets_[fieldIndex]; + const auto& worksetLIDs = worksetLIDs_; + const auto& fieldValues = gatherFields_[fieldIndex]; + + Kokkos::parallel_for(Kokkos::RangePolicy(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) { + for(int basis=0; basis < static_cast(fieldOffsets.size()); ++basis) { + const int lid = worksetLIDs(cell,fieldOffsets(basis)); + fieldValues(cell,basis) = kokkosSolution(lid,0); } + }); + } - // loop over the fields to be gathered - Teuchos::ArrayRCP local_x; - for (std::size_t fieldIndex=0; fieldIndexgetFieldBlock(fieldNum); - - // grab local data for inputing - RCP > block_x = rcp_dynamic_cast >(x->getNonconstVectorBlock(indexerId)); - block_x->getLocalData(ptrFromRef(local_x)); - - const std::vector & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum); - - // loop over basis functions and fill the fields - for(std::size_t basis=0;basis ScatterResidual_BlockedTpetra(); }; +// ************************************************************** +// Tangent +// ************************************************************** +template +class ScatterResidual_BlockedTpetra + : public panzer::EvaluatorWithBaseImpl, + public PHX::EvaluatorDerived, + public panzer::CloneableEvaluator { + +public: + ScatterResidual_BlockedTpetra(const Teuchos::RCP & indexer) + : globalIndexer_(indexer) {} + + ScatterResidual_BlockedTpetra(const Teuchos::RCP & indexer, + const Teuchos::ParameterList& p); + + void postRegistrationSetup(typename TRAITS::SetupData d, + PHX::FieldManager& vm); + + void preEvaluate(typename TRAITS::PreEvalData d); + + void evaluateFields(typename TRAITS::EvalData workset); + + virtual Teuchos::RCP clone(const Teuchos::ParameterList & pl) const + { return Teuchos::rcp(new ScatterResidual_BlockedTpetra(globalIndexer_,pl)); } + +private: + typedef typename panzer::Traits::Tangent::ScalarT ScalarT; + typedef typename TRAITS::RealType RealType; + + typedef BlockedTpetraLinearObjContainer ContainerType; + typedef Tpetra::Vector VectorType; + typedef Tpetra::CrsMatrix CrsMatrixType; + typedef Tpetra::CrsGraph CrsGraphType; + typedef Tpetra::Map MapType; + typedef Tpetra::Import ImportType; + typedef Tpetra::Export ExportType; + + //! Dummy evalauted field so that the evaluator will have something to do + Teuchos::RCP scatterHolder_; + + //! Fields that need to be scattered will be put in this vector + std::vector< PHX::MDField > scatterFields_; + + //! Maps the local (field,element,basis) triplet to a global ID for scattering + Teuchos::RCP globalIndexer_; + + //! Vector of global indexers, one for each scattered field respectively + std::vector> fieldGlobalIndexers_; + + //! Field IDs in the local product vector block (not global field id) + std::vector fieldIds_; + + //! Returns the index into the Thyra ProductVector sub-block. Size + //! of number of fields to scatter + std::vector productVectorBlockIndex_; + + // This maps the scattered field names to the DOF manager field + // For instance a Navier-Stokes map might look like + // fieldMap_["RESIDUAL_Velocity"] --> "Velocity" + // fieldMap_["RESIDUAL_Pressure"] --> "Pressure" + Teuchos::RCP > fieldMap_; + + std::string globalDataKey_; // what global data does this fill? + Teuchos::RCP > blockedContainer_; + + /// Storage for the tangent data + PHX::ViewOfViews<2,Kokkos::View> dfdpFieldsVoV_; + + //! Local indices for unknowns + PHX::View worksetLIDs_; + + //! Offset into the cell lids for each field + std::vector> fieldOffsets_; + + ScatterResidual_BlockedTpetra(); +}; + } #ifdef Panzer_BUILD_HESSIAN_SUPPORT diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedTpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedTpetra_impl.hpp index 743d8b30cef7..3bb426ec0d34 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedTpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_BlockedTpetra_impl.hpp @@ -22,6 +22,7 @@ #include "Panzer_BlockedTpetraLinearObjContainer.hpp" #include "Panzer_LOCPair_GlobalEvaluationData.hpp" #include "Panzer_HashUtils.hpp" +#include "Panzer_ParameterList_GlobalEvaluationData.hpp" #include "Panzer_GlobalEvaluationDataContainer.hpp" #include "Thyra_ProductVectorBase.hpp" @@ -472,4 +473,171 @@ evaluateFields(typename TRAITS::EvalData workset) } +// ********************************************************************** +// Specialization: Tangent +// ********************************************************************** + +template +panzer::ScatterResidual_BlockedTpetra:: +ScatterResidual_BlockedTpetra(const Teuchos::RCP & indexer, + const Teuchos::ParameterList& p) + : globalIndexer_(indexer) + , globalDataKey_("Residual Scatter Container") +{ + std::string scatterName = p.get("Scatter Name"); + scatterHolder_ = + Teuchos::rcp(new PHX::Tag(scatterName,Teuchos::rcp(new PHX::MDALayout(0)))); + + // get names to be evaluated + const std::vector& names = + *(p.get< Teuchos::RCP< std::vector > >("Dependent Names")); + + // grab map from evaluated names to field names + fieldMap_ = p.get< Teuchos::RCP< std::map > >("Dependent Map"); + + Teuchos::RCP dl = + p.get< Teuchos::RCP >("Basis")->functional; + + // build the vector of fields that this is dependent on + scatterFields_.resize(names.size()); + for (std::size_t eq = 0; eq < names.size(); ++eq) { + scatterFields_[eq] = PHX::MDField(names[eq],dl); + + // tell the field manager that we depend on this field + this->addDependentField(scatterFields_[eq]); + } + + // this is what this evaluator provides + this->addEvaluatedField(*scatterHolder_); + + if (p.isType("Global Data Key")) + globalDataKey_ = p.get("Global Data Key"); + + this->setName(scatterName+" Scatter Residual (Tangent)"); +} + +// ********************************************************************** +template +void panzer::ScatterResidual_BlockedTpetra:: +postRegistrationSetup(typename TRAITS::SetupData d, + PHX::FieldManager& /* fm */) +{ + const Workset & workset_0 = (*d.worksets_)[0]; + const std::string blockId = this->wda(workset_0).block_id; + + fieldIds_.resize(scatterFields_.size()); + fieldOffsets_.resize(scatterFields_.size()); + fieldGlobalIndexers_.resize(scatterFields_.size()); + productVectorBlockIndex_.resize(scatterFields_.size()); + int maxElementBlockGIDCount = -1; + for(std::size_t fd=0; fd < scatterFields_.size(); ++fd) { + const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second; + const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager + productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum); + fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]]; + fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName); // Field number in the sub-global-indexer + + const std::vector& offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]); + fieldOffsets_[fd] = PHX::View("ScatterResidual_BlockedTpetra(Tangent):fieldOffsets",offsets.size()); + auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]); + for (std::size_t i=0; i < offsets.size(); ++i) + hostOffsets(i) = offsets[i]; + Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets); + + maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount); + } + + // We will use one workset lid view for all fields, but has to be + // sized big enough to hold the largest elementBlockGIDCount in the + // ProductVector. + worksetLIDs_ = PHX::View("ScatterResidual_BlockedTpetra(Tangent):worksetLIDs", + scatterFields_[0].extent(0), + maxElementBlockGIDCount); +} + +// ********************************************************************** +template +void panzer::ScatterResidual_BlockedTpetra:: +preEvaluate(typename TRAITS::PreEvalData d) +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Thyra::ProductVectorBase; + + // this is the list of parameters and their names that this scatter has to account for + std::vector activeParameters = + rcp_dynamic_cast(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters(); + + const int numBlocks = static_cast(globalIndexer_->getFieldDOFManagers().size()); + + dfdpFieldsVoV_.initialize("ScatterResidual_Tpetra::dfdpFieldsVoV_",activeParameters.size(),numBlocks); + + for(std::size_t i=0;i paramBlockedContainer = rcp_dynamic_cast(d.gedc->getDataObject(activeParameters[i]),true); + RCP> productVector = + rcp_dynamic_cast>(paramBlockedContainer->get_f(),true); + for(int j=0;j>(productVector->getNonconstVectorBlock(j),true))->getTpetraVector()); + const auto& dfdp_view = tpetraBlock.getLocalViewDevice(Tpetra::Access::ReadWrite); + dfdpFieldsVoV_.addView(dfdp_view,i,j); + } + } + + dfdpFieldsVoV_.syncHostToDevice(); + + // extract linear object container + blockedContainer_ = rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_)); + + if(blockedContainer_==Teuchos::null) { + RCP gdata = rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_),true); + blockedContainer_ = rcp_dynamic_cast(gdata->getGhostedLOC()); + } +} + +// ********************************************************************** +template +void panzer::ScatterResidual_BlockedTpetra:: +evaluateFields(typename TRAITS::EvalData workset) +{ + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using Thyra::VectorBase; + using Thyra::ProductVectorBase; + + const auto& localCellIds = this->wda(workset).cell_local_ids_k; + const RCP> thyraBlockResidual = rcp_dynamic_cast >(blockedContainer_->get_f(),true); + + // Loop over scattered fields + int currentWorksetLIDSubBlock = -1; + for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { + // workset LIDs only change for different sub blocks + if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) { + fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_); + currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex]; + } + + auto& tpetraResidual = *((rcp_dynamic_cast>(thyraBlockResidual->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),true))->getTpetraVector()); + const auto& kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::ReadWrite); + + // Class data fields for lambda capture + const auto& fieldOffsets = fieldOffsets_[fieldIndex]; + const auto& worksetLIDs = worksetLIDs_; + const auto& fieldValues = scatterFields_[fieldIndex].get_static_view(); + const auto& tangentFieldsDevice = dfdpFieldsVoV_.getViewDevice(); + const auto& kokkosTangents = Kokkos::subview(tangentFieldsDevice,Kokkos::ALL(),productVectorBlockIndex_[fieldIndex]); + const double num_params = Kokkos::dimension_scalar(fieldValues)-1; + + Kokkos::parallel_for(Kokkos::RangePolicy(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) { + for(int basis=0; basis < static_cast(fieldOffsets.size()); ++basis) { + const int lid = worksetLIDs(cell,fieldOffsets(basis)); + Kokkos::atomic_add(&kokkosResidual(lid,0), fieldValues(cell,basis).val()); + for(int i_param=0; i_param