From 4bb2336e5bb442d41b83e0077b60c86e3e04c00d Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Wed, 6 Nov 2024 12:59:06 -0700 Subject: [PATCH 01/11] Panzer Tangents :: Small progress for scatter, move unit test Signed-off-by: Bryan Reuter --- .../test/evaluator_tests/CMakeLists.txt | 7 + .../tpetra_blocked_scatter_residual.cpp | 546 ++++++++++++++++++ .../tpetra_scatter_residual.cpp | 136 ++--- .../Panzer_ScatterResidual_Tpetra_impl.hpp | 100 +++- 4 files changed, 669 insertions(+), 120 deletions(-) create mode 100644 packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt b/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt index d871d0375cb0..079cdad0848b 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt +++ b/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt @@ -57,6 +57,13 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( NUM_MPI_PROCS 2 ) +TRIBITS_ADD_EXECUTABLE_AND_TEST( + tScatterResidual_BlockedTpetra + SOURCES tpetra_blocked_scatter_residual.cpp ${UNIT_TEST_DRIVER} + COMM serial mpi + NUM_MPI_PROCS 2 + ) + TRIBITS_ADD_EXECUTABLE_AND_TEST( tScatterDirichletResidual_Tpetra SOURCES tpetra_scatter_dirichlet_residual.cpp ${UNIT_TEST_DRIVER} 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 new file mode 100644 index 000000000000..a5e3d75d2362 --- /dev/null +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp @@ -0,0 +1,546 @@ +// @HEADER +// ***************************************************************************** +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// +// Copyright 2011 NTESS and the Panzer contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#include +#include +#include +#include + +using Teuchos::RCP; +using Teuchos::rcp; + +#include "Teuchos_DefaultComm.hpp" +#include "Teuchos_GlobalMPISession.hpp" + +#include "Panzer_FieldManagerBuilder.hpp" +#include "Panzer_BlockedDOFManager.hpp" +#include "Panzer_BlockedTpetraLinearObjFactory.hpp" +#include "Panzer_PureBasis.hpp" +#include "Panzer_BasisIRLayout.hpp" +#include "Panzer_Workset.hpp" +#include "Panzer_GatherOrientation.hpp" +#include "Panzer_ScatterResidual_BlockedTpetra.hpp" +#include "Panzer_GatherSolution_BlockedTpetra.hpp" +#include "Panzer_GlobalEvaluationDataContainer.hpp" + +#include "Panzer_STK_Version.hpp" +#include "PanzerAdaptersSTK_config.hpp" +#include "Panzer_STK_Interface.hpp" +#include "Panzer_STK_SquareQuadMeshFactory.hpp" +#include "Panzer_STK_SetupUtilities.hpp" +#include "Panzer_STKConnManager.hpp" + +#include "Teuchos_DefaultMpiComm.hpp" +#include "Teuchos_OpaqueWrapper.hpp" + +#include "Thyra_VectorStdOps.hpp" +#include "Thyra_ProductVectorBase.hpp" +#include "Thyra_SpmdVectorBase.hpp" + +#include "user_app_EquationSetFactory.hpp" + +#include // for get char +#include +#include + +#include "Panzer_Evaluator_WithBaseImpl.hpp" + +namespace panzer +{ + typedef Teuchos::ArrayRCP::size_type size_type; + + using TpetraBlockedLinObjFactoryType = panzer::BlockedTpetraLinearObjFactory; + using TpetraBlockedLinObjContainerType = panzer::BlockedTpetraLinearObjContainer; + + Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName); + void testInitialization(const Teuchos::RCP &ipb); + Teuchos::RCP buildMesh(int elemX, int elemY); + + TEUCHOS_UNIT_TEST(block_assembly, scatter_solution_residual) + { + +#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"; + + 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); + + // setup field manager, add evaluator under test + ///////////////////////////////////////////////////////////// + + PHX::FieldManager fm; + + 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> evaluator = lof->buildGather(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> evaluator = lof->buildGather(pl); + + fm.registerEvaluator(evaluator); + } + + std::vector derivative_dimensions; + derivative_dimensions.push_back(12); + fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); + + 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); + 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); + + // test Residual 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); + } + } + + TEUCHOS_UNIT_TEST(block_assembly, scatter_solution_jacobian) + { + +#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"; + + 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 | LinearObjContainer::Mat, *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); + + // setup field manager, add evaluator under test + ///////////////////////////////////////////////////////////// + + PHX::FieldManager fm; + + 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> evaluator = lof->buildGather(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> evaluator = lof->buildGather(pl); + + fm.registerEvaluator(evaluator); + } + + std::vector derivative_dimensions; + derivative_dimensions.push_back(12); + fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); + + 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); + 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); + + // test Jacobian 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()); + out << std::endl; + for (size_type i = 0; i < data.size(); i++) + { + double target = 123.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + out << data[i] << std::endl; + } + + Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(1))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(1)->getLocalNumElements()); + out << std::endl; + for (size_type i = 0; i < data.size(); i++) + { + double target = 456.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + out << data[i] << std::endl; + } + + Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(2))->getLocalData(Teuchos::ptrFromRef(data)); + TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(2)->getLocalNumElements()); + out << std::endl; + for (size_type i = 0; i < data.size(); i++) + { + double target = 789.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + out << data[i] << std::endl; + } + } + + Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName) + { + Teuchos::RCP topo = + Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData>())); + + panzer::CellData cellData(worksetSize, topo); + return Teuchos::rcp(new panzer::PureBasis(basisName, 1, cellData)); + } + + Teuchos::RCP buildMesh(int elemX, int elemY) + { + RCP pl = rcp(new Teuchos::ParameterList); + pl->set("X Blocks", 1); + pl->set("Y Blocks", 1); + pl->set("X Elements", elemX); + pl->set("Y Elements", elemY); + + panzer_stk::SquareQuadMeshFactory factory; + factory.setParameterList(pl); + RCP mesh = factory.buildUncommitedMesh(MPI_COMM_WORLD); + factory.completeMeshConstruction(*mesh, MPI_COMM_WORLD); + + return mesh; + } + + void testInitialization(const Teuchos::RCP &ipb) + { + // Physics block + ipb->setName("test physics"); + { + Teuchos::ParameterList &p = ipb->sublist("a"); + p.set("Type", "Energy"); + p.set("Prefix", ""); + p.set("Model ID", "solid"); + p.set("Basis Type", "HGrad"); + p.set("Basis Order", 1); + p.set("Integration Order", 1); + } + { + Teuchos::ParameterList &p = ipb->sublist("b"); + p.set("Type", "Energy"); + p.set("Prefix", "ION_"); + p.set("Model ID", "solid"); + p.set("Basis Type", "HCurl"); + p.set("Basis Order", 1); + p.set("Integration Order", 1); + } + } + +} diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp index a5e3d75d2362..27d03219961b 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp @@ -20,14 +20,14 @@ using Teuchos::rcp; #include "Teuchos_GlobalMPISession.hpp" #include "Panzer_FieldManagerBuilder.hpp" -#include "Panzer_BlockedDOFManager.hpp" -#include "Panzer_BlockedTpetraLinearObjFactory.hpp" +#include "Panzer_DOFManager.hpp" +#include "Panzer_TpetraLinearObjFactory.hpp" #include "Panzer_PureBasis.hpp" #include "Panzer_BasisIRLayout.hpp" #include "Panzer_Workset.hpp" #include "Panzer_GatherOrientation.hpp" -#include "Panzer_ScatterResidual_BlockedTpetra.hpp" -#include "Panzer_GatherSolution_BlockedTpetra.hpp" +#include "Panzer_ScatterResidual_Tpetra.hpp" +#include "Panzer_GatherSolution_Tpetra.hpp" #include "Panzer_GlobalEvaluationDataContainer.hpp" #include "Panzer_STK_Version.hpp" @@ -41,7 +41,7 @@ using Teuchos::rcp; #include "Teuchos_OpaqueWrapper.hpp" #include "Thyra_VectorStdOps.hpp" -#include "Thyra_ProductVectorBase.hpp" +#include "Thyra_VectorBase.hpp" #include "Thyra_SpmdVectorBase.hpp" #include "user_app_EquationSetFactory.hpp" @@ -56,14 +56,16 @@ namespace panzer { typedef Teuchos::ArrayRCP::size_type size_type; - using TpetraBlockedLinObjFactoryType = panzer::BlockedTpetraLinearObjFactory; - using TpetraBlockedLinObjContainerType = panzer::BlockedTpetraLinearObjContainer; + using TpetraLinObjFactoryType = panzer::TpetraLinearObjFactory; + using TpetraLinObjContainerType = panzer::TpetraLinearObjContainer; Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName); void testInitialization(const Teuchos::RCP &ipb); Teuchos::RCP buildMesh(int elemX, int elemY); - TEUCHOS_UNIT_TEST(block_assembly, scatter_solution_residual) + // TODO BWR NO TANGENT COVERAGE HERE OR IN OTHER TEST + + TEUCHOS_UNIT_TEST(assembly, scatter_solution_residual) { #ifdef HAVE_MPI @@ -102,35 +104,32 @@ namespace panzer // 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)); + RCP dofManager = Teuchos::rcp(new panzer::DOFManager(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); + std::vector fieldOrder; + fieldOrder.push_back(fieldName1_q1); + fieldOrder.push_back(fieldName_qedge1); + fieldOrder.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); + Teuchos::RCP t_lof = Teuchos::rcp(new TpetraLinObjFactoryType(tComm.getConst(), dofManager)); + Teuchos::RCP> lof = t_lof; + Teuchos::RCP loc = t_lof->buildGhostedLinearObjContainer(); + t_lof->initializeGhostedContainer(LinearObjContainer::X | LinearObjContainer::F, *loc); loc->initialize(); - Teuchos::RCP b_loc = Teuchos::rcp_dynamic_cast(loc); + Teuchos::RCP t_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); + Teuchos::RCP> x_vec = t_loc->get_x_th(); + Thyra::assign(x_vec.ptr(), 123.0 + myRank); // setup field manager, add evaluator under test ///////////////////////////////////////////////////////////// @@ -218,10 +217,6 @@ namespace panzer fm.registerEvaluator(evaluator); } - std::vector derivative_dimensions; - derivative_dimensions.push_back(12); - fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); - panzer::Traits::SD sd; sd.worksets_ = work_sets; @@ -242,41 +237,27 @@ namespace panzer workset.evaluate_transient_terms = false; fm.evaluateFields(workset); + fm.postEvaluate(0); // test Residual fields Teuchos::ArrayRCP data; - Teuchos::RCP> f_vec = Teuchos::rcp_dynamic_cast>(b_loc->get_f()); + Teuchos::RCP> f_vec = t_loc->get_f_th(); // 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()); + Teuchos::rcp_dynamic_cast>(f_vec)->getLocalData(Teuchos::ptrFromRef(data)); 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); - } } - TEUCHOS_UNIT_TEST(block_assembly, scatter_solution_jacobian) + TEUCHOS_UNIT_TEST(assembly, scatter_solution_jacobian) { + // TODO BWR not checking the jacobian calculation, just the residual part!!! #ifdef HAVE_MPI Teuchos::RCP> tComm = Teuchos::rcp(new Teuchos::MpiComm(MPI_COMM_WORLD)); @@ -314,42 +295,42 @@ namespace panzer // 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)); + RCP dofManager = Teuchos::rcp(new panzer::DOFManager(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); + std::vector fieldOrder; + fieldOrder.push_back(fieldName1_q1); + fieldOrder.push_back(fieldName_qedge1); + fieldOrder.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 | LinearObjContainer::Mat, *loc); + Teuchos::RCP t_lof = Teuchos::rcp(new TpetraLinObjFactoryType(tComm.getConst(), dofManager)); + Teuchos::RCP> lof = t_lof; + Teuchos::RCP loc = t_lof->buildGhostedLinearObjContainer(); + t_lof->initializeGhostedContainer(LinearObjContainer::X | LinearObjContainer::F | LinearObjContainer::Mat, *loc); loc->initialize(); - Teuchos::RCP b_loc = Teuchos::rcp_dynamic_cast(loc); + Teuchos::RCP t_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); + Teuchos::RCP> x_vec = t_loc->get_x_th(); + Thyra::assign(x_vec.ptr(), 123.0 + myRank); // setup field manager, add evaluator under test ///////////////////////////////////////////////////////////// PHX::FieldManager fm; + std::vector derivative_dimensions; + derivative_dimensions.push_back(12); + fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); + std::string resName = ""; Teuchos::RCP> names_map = Teuchos::rcp(new std::map); @@ -431,10 +412,6 @@ namespace panzer fm.registerEvaluator(evaluator); } - std::vector derivative_dimensions; - derivative_dimensions.push_back(12); - fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); - panzer::Traits::SD sd; sd.worksets_ = work_sets; @@ -455,43 +432,22 @@ namespace panzer workset.evaluate_transient_terms = false; fm.evaluateFields(workset); + fm.postEvaluate(0); // test Jacobian fields Teuchos::ArrayRCP data; - Teuchos::RCP> f_vec = Teuchos::rcp_dynamic_cast>(b_loc->get_f()); + Teuchos::RCP> f_vec = t_loc->get_f_th(); // 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()); - out << std::endl; + Teuchos::rcp_dynamic_cast>(f_vec)->getLocalData(Teuchos::ptrFromRef(data)); for (size_type i = 0; i < data.size(); i++) { double target = 123.0 + myRank; TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); - out << data[i] << std::endl; - } - - Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(1))->getLocalData(Teuchos::ptrFromRef(data)); - TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(1)->getLocalNumElements()); - out << std::endl; - for (size_type i = 0; i < data.size(); i++) - { - double target = 456.0 + myRank; - TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); - out << data[i] << std::endl; } - Teuchos::rcp_dynamic_cast>(f_vec->getVectorBlock(2))->getLocalData(Teuchos::ptrFromRef(data)); - TEST_EQUALITY(static_cast(data.size()), b_loc->getMapForBlock(2)->getLocalNumElements()); - out << std::endl; - for (size_type i = 0; i < data.size(); i++) - { - double target = 789.0 + myRank; - TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); - out << data[i] << std::endl; - } } Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName) diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 117f271faa56..0a197e01fe17 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -146,6 +146,7 @@ ScatterResidual_Tpetra(const Teuchos::RCP & indexer // build the vector of fields that this is dependent on scatterFields_.resize(names.size()); + scratch_offsets_.resize(names.size()); for (std::size_t eq = 0; eq < names.size(); ++eq) { scatterFields_[eq] = PHX::MDField(names[eq],dl); @@ -169,12 +170,21 @@ postRegistrationSetup(typename TRAITS::SetupData /* d */, PHX::FieldManager& /* fm */) { fieldIds_.resize(scatterFields_.size()); + const Workset & workset_0 = (*d.worksets_)[0]; + std::string blockId = this->wda(workset_0).block_id; + // load required field numbers for fast use for(std::size_t fd=0;fdfind(scatterFields_[fd].fieldTag().name())->second; fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName); + + const std::vector & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]); + scratch_offsets_[fd] = PHX::View("offsets",offsets.size()); + Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View(offsets.data(), offsets.size())); } + scratch_lids_ = PHX::View("lids",scatterFields_[0].extent(0), + globalIndexer_->getElementBlockGIDCount(blockId)); } // ********************************************************************** @@ -195,9 +205,21 @@ preEvaluate(typename TRAITS::PreEvalData d) for(std::size_t i=0;i vec = rcp_dynamic_cast(d.gedc->getDataObject(activeParameters[i]),true)->get_f(); + auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access:ReadOnly); + Teuchos::ArrayRCP vec_array = vec->get1dViewNonConst(); dfdp_vectors_.push_back(vec_array); } + + // TODO BWR DO WE WANT TO BE ABLE TO FILL RESIDUAL TOO? + // extract linear object container + tpetraContainer_ = Teuchos::rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_)); + + if(tpetraContainer_==Teuchos::null) { + // extract linear object container + Teuchos::RCP loc = Teuchos::rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC(); + tpetraContainer_ = Teuchos::rcp_dynamic_cast(loc); + } } // ********************************************************************** @@ -205,36 +227,54 @@ template void panzer::ScatterResidual_Tpetra:: evaluateFields(typename TRAITS::EvalData workset) { - // 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; - - // NOTE: A reordering of these loops will likely improve performance - // The "getGIDFieldOffsets may be expensive. However the - // "getElementGIDs" can be cheaper. However the lookup for LIDs - // may be more expensive! - - // scatter operation for each cell in workset - for(std::size_t worksetCellIndex=0;worksetCellIndexgetElementLIDs(cellLocalId); - - // loop over each field to be scattered - for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { - int fieldNum = fieldIds_[fieldIndex]; - const std::vector & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum); - - // loop over basis functions - for(std::size_t basis=0;basis LOC; + + // for convenience pull out some objects from workset + std::string blockId = this->wda(workset).block_id; + + Teuchos::RCP r = tpetraContainer_->get_f(); + + globalIndexer_->getElementLIDs(this->wda(workset).getLocalCellIDs(),scratch_lids_); + + ScatterResidual_Tangent_Functor functor; + functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite); + functor.lids = scratch_lids_; + + // for each field, do a parallel for loop + for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { + functor.offsets = scratch_offsets_[fieldIndex]; + functor.field = scatterFields_[fieldIndex]; + functor.tangents = dfdp_vectors_; // TODO ??? + + Kokkos::parallel_for(workset.num_cells,functor); + } + + // NOTE: A reordering of these loops will likely improve performance + // The "getGIDFieldOffsets may be expensive. However the + // "getElementGIDs" can be cheaper. However the lookup for LIDs + // may be more expensive! + + // scatter operation for each cell in workset + for(std::size_t worksetCellIndex=0;worksetCellIndexgetElementLIDs(cellLocalId); + + // loop over each field to be scattered + for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { + int fieldNum = fieldIds_[fieldIndex]; + const std::vector & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum); + + // loop over basis functions + for(std::size_t basis=0;basis Date: Wed, 6 Nov 2024 15:59:55 -0700 Subject: [PATCH 02/11] Panzer Tangents :: First pass at implementation Signed-off-by: Bryan Reuter --- ...nzer_GatherSolution_BlockedTpetra_impl.hpp | 3 - .../Panzer_ScatterResidual_Tpetra_decl.hpp | 9 +- .../Panzer_ScatterResidual_Tpetra_impl.hpp | 141 ++++++++++-------- 3 files changed, 87 insertions(+), 66 deletions(-) diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp index 52488585d37e..9c8533419013 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp @@ -305,9 +305,6 @@ postRegistrationSetup(typename TRAITS::SetupData d, maxElementBlockGIDCount); // Set up storage for tangentFields using view of views - // We also need storage for the number of tangent fields associated with - // each gatherField - if (has_tangent_fields_) { size_t inner_vector_max_size = 0; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp index 4ea8c009ed5c..04609a23a076 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp @@ -14,6 +14,7 @@ #include "Phalanx_config.hpp" #include "Phalanx_Evaluator_Macros.hpp" #include "Phalanx_MDField.hpp" +#include "Phalanx_KokkosViewOfViews.hpp" #include "Teuchos_ParameterList.hpp" @@ -134,6 +135,7 @@ class ScatterResidual_Tpetra private: typedef typename panzer::Traits::Tangent::ScalarT ScalarT; + typedef typename panzer::Traits::RealType RealT; // dummy field so that the evaluator will have something to do Teuchos::RCP scatterHolder_; @@ -153,8 +155,13 @@ class ScatterResidual_Tpetra Teuchos::RCP > fieldMap_; std::string globalDataKey_; // what global data does this fill? + Teuchos::RCP > tpetraContainer_; + + /// Storage for the tangent data + PHX::ViewOfViews<1,PHX::View> dfdpFieldsVoV_; - std::vector< Teuchos::ArrayRCP > dfdp_vectors_; + PHX::View scratch_lids_; + std::vector > scratch_offsets_; ScatterResidual_Tpetra(); }; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 0a197e01fe17..7639dc582391 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -24,6 +24,7 @@ #include "Panzer_GlobalEvaluationDataContainer.hpp" #include "Phalanx_DataLayout_MDALayout.hpp" +#include "Phalanx_Scratch_Utilities.hpp" #include "Teuchos_FancyOStream.hpp" @@ -130,6 +131,7 @@ ScatterResidual_Tpetra(const Teuchos::RCP & indexer : 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)))); @@ -166,7 +168,7 @@ ScatterResidual_Tpetra(const Teuchos::RCP & indexer // ********************************************************************** template void panzer::ScatterResidual_Tpetra:: -postRegistrationSetup(typename TRAITS::SetupData /* d */, +postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager& /* fm */) { fieldIds_.resize(scatterFields_.size()); @@ -201,16 +203,20 @@ preEvaluate(typename TRAITS::PreEvalData d) std::vector activeParameters = rcp_dynamic_cast(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters(); - dfdp_vectors_.clear(); + dfdpFieldsVoV_.initialize("ScatterResidual_Tpetra::dfdpFieldsVoV_",activeParameters.size()); + for(std::size_t i=0;i vec = rcp_dynamic_cast(d.gedc->getDataObject(activeParameters[i]),true)->get_f(); - auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access:ReadOnly); - - Teuchos::ArrayRCP vec_array = vec->get1dViewNonConst(); - dfdp_vectors_.push_back(vec_array); + auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access::ReadWrite); + + // TODO BWR not sure this will be checked by tests... + TEUCHOS_ASSERT(dfdp_view.extent_int(1)==1); + dfdpFieldsVoV_.addView(Kokkos::subview(dfdp_view,Kokkos::ALL(),0),i); } + dfdpFieldsVoV_.syncHostToDevice(); + // TODO BWR DO WE WANT TO BE ABLE TO FILL RESIDUAL TOO? // extract linear object container tpetraContainer_ = Teuchos::rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_)); @@ -222,61 +228,6 @@ preEvaluate(typename TRAITS::PreEvalData d) } } -// ********************************************************************** -template -void panzer::ScatterResidual_Tpetra:: -evaluateFields(typename TRAITS::EvalData workset) -{ - typedef TpetraLinearObjContainer LOC; - - // for convenience pull out some objects from workset - std::string blockId = this->wda(workset).block_id; - - Teuchos::RCP r = tpetraContainer_->get_f(); - - globalIndexer_->getElementLIDs(this->wda(workset).getLocalCellIDs(),scratch_lids_); - - ScatterResidual_Tangent_Functor functor; - functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite); - functor.lids = scratch_lids_; - - // for each field, do a parallel for loop - for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { - functor.offsets = scratch_offsets_[fieldIndex]; - functor.field = scatterFields_[fieldIndex]; - functor.tangents = dfdp_vectors_; // TODO ??? - - Kokkos::parallel_for(workset.num_cells,functor); - } - - // NOTE: A reordering of these loops will likely improve performance - // The "getGIDFieldOffsets may be expensive. However the - // "getElementGIDs" can be cheaper. However the lookup for LIDs - // may be more expensive! - - // scatter operation for each cell in workset - for(std::size_t worksetCellIndex=0;worksetCellIndexgetElementLIDs(cellLocalId); - - // loop over each field to be scattered - for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { - int fieldNum = fieldIds_[fieldIndex]; - const std::vector & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum); - - // loop over basis functions - for(std::size_t basis=0;basis offsets; // how to get a particular field FieldType field; - KOKKOS_INLINE_FUNCTION void operator()(const unsigned int cell) const { @@ -447,6 +397,45 @@ class ScatterResidual_Residual_Functor { } }; +template +class ScatterResidual_Tangent_Functor { +public: + typedef typename PHX::Device execution_space; + typedef PHX::MDField FieldType; + + bool fillResidual; + Kokkos::View r_data; + + Kokkos::View lids; // local indices for unknowns. + PHX::View offsets; // how to get a particular field + FieldType field; + + PHX::ViewOfViews<1,PHX::View> dfdp_fields; // tangent fields + + KOKKOS_INLINE_FUNCTION + void operator()(const unsigned int cell) const + { + + const std::size_t numDerivs = PHX::getFadSize(field.get_static_view()); + const auto tangents = dfdp_fields.getViewDevice(); + // loop over the basis functions (currently they are nodes) + for(std::size_t basis=0; basis < offsets.extent(0); basis++) { + typename FieldType::array_type::reference_type scatterField = field(cell,basis); + int offset = offsets(basis); + LO lid = lids(cell,offset); + + // Sum residual + if(fillResidual) + Kokkos::atomic_add(&r_data(lid,0), scatterField.val()); + + // loop over the tangents + for(int i_param=0; i_param +void panzer::ScatterResidual_Tpetra:: +evaluateFields(typename TRAITS::EvalData workset) +{ + typedef TpetraLinearObjContainer LOC; + + // for convenience pull out some objects from workset + std::string blockId = this->wda(workset).block_id; + + Teuchos::RCP r = tpetraContainer_->get_f(); + + globalIndexer_->getElementLIDs(this->wda(workset).getLocalCellIDs(),scratch_lids_); + + ScatterResidual_Tangent_Functor functor; + functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite); + functor.lids = scratch_lids_; + + // for each field, do a parallel for loop + for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { + functor.offsets = scratch_offsets_[fieldIndex]; + functor.field = scatterFields_[fieldIndex]; + functor.dfdp_fields = dfdpFieldsVoV_; + + Kokkos::parallel_for(workset.num_cells,functor); + } +} + // ********************************************************************** #endif From 8d2a055910582c8f554df9cf5b1fe8ad78adb210 Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Thu, 7 Nov 2024 17:55:08 -0700 Subject: [PATCH 03/11] Panzer Tangents :: On the right track but my test is not quite set up correctly Signed-off-by: Bryan Reuter --- .../tpetra_scatter_residual.cpp | 314 ++++++++++++++++++ .../Panzer_ScatterResidual_Tpetra_impl.hpp | 17 +- 2 files changed, 324 insertions(+), 7 deletions(-) diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp index 27d03219961b..3c63810a0cd5 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp @@ -29,6 +29,8 @@ using Teuchos::rcp; #include "Panzer_ScatterResidual_Tpetra.hpp" #include "Panzer_GatherSolution_Tpetra.hpp" #include "Panzer_GlobalEvaluationDataContainer.hpp" +#include "Panzer_LOCPair_GlobalEvaluationData.hpp" +#include "Panzer_ParameterList_GlobalEvaluationData.hpp" #include "Panzer_STK_Version.hpp" #include "PanzerAdaptersSTK_config.hpp" @@ -450,6 +452,318 @@ namespace panzer } + TEUCHOS_UNIT_TEST(assembly, scatter_solution_tangent) + { + // TODO BWR not checking the tangent calculation, just the residual part!!! + +#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"; + + 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::DOFManager(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; + fieldOrder.push_back(fieldName1_q1); + fieldOrder.push_back(fieldName_qedge1); + fieldOrder.push_back(fieldName2_q1); + dofManager->setFieldOrder(fieldOrder); + + dofManager->buildGlobalUnknowns(); + + // setup linear object factory + ///////////////////////////////////////////////////////////// + Teuchos::RCP t_lof = Teuchos::rcp(new TpetraLinObjFactoryType(tComm.getConst(), dofManager)); + Teuchos::RCP> lof = t_lof; + Teuchos::RCP loc = t_lof->buildGhostedLinearObjContainer(); + t_lof->initializeGhostedContainer(LinearObjContainer::X | LinearObjContainer::F , *loc); + loc->initialize(); + + Teuchos::RCP t_loc = Teuchos::rcp_dynamic_cast(loc); + Teuchos::RCP> x_vec = t_loc->get_x_th(); + Thyra::assign(x_vec.ptr(), 123.0 + myRank); + + std::vector> tangentContainers; + + using LOCPair = panzer::LOCPair_GlobalEvaluationData; + using Teuchos::rcp_dynamic_cast; + + auto locPair = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X)); + + auto global_t_loc = rcp_dynamic_cast(locPair->getGlobalLOC()); + Teuchos::RCP> global_x_vec = global_t_loc->get_x_th(); + Thyra::assign(global_x_vec.ptr(), 0.123 + myRank); + + auto ghosted_t_loc = rcp_dynamic_cast(locPair->getGhostedLOC()); + Teuchos::RCP> ghosted_x_vec = ghosted_t_loc->get_x_th(); + Thyra::assign(ghosted_x_vec.ptr(), 0.123 + myRank); + + tangentContainers.push_back(locPair); + + // setup field manager, add evaluator under test + ///////////////////////////////////////////////////////////// + + PHX::FieldManager fm; + + std::vector derivative_dimensions; + derivative_dimensions.push_back(1); + 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)); + (*tangent_names)[0].push_back(fieldName1_q1 + " Tangent"); + (*tangent_names)[1].push_back(fieldName2_q1 + " Tangent"); + pl.set("Tangent Names", tangent_names); + + Teuchos::RCP> evaluator = lof->buildGather(pl); + + fm.registerEvaluator(evaluator); + } + { + 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); + tangent_names->push_back(fieldName1_q1 + " Tangent"); + tangent_names->push_back(fieldName2_q1 + " Tangent"); + + Teuchos::ParameterList pl; + pl.set("Basis", basis_q1); + pl.set("DOF Names", tangent_names); + pl.set("Indexer Names", names); + + pl.set("Global Data Key", "Tangent Container"); + + 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)); + (*tangent_names)[0].push_back(fieldName_qedge1 + " Tangent"); + pl.set("Tangent Names", tangent_names); + + Teuchos::RCP> evaluator = lof->buildGather(pl); + + fm.registerEvaluator(evaluator); + } + { + using Teuchos::RCP; + using Teuchos::rcp; + RCP> names = rcp(new std::vector); + RCP> tangent_names = rcp(new std::vector); + names->push_back(fieldName_qedge1); + tangent_names->push_back(fieldName_qedge1 + " Tangent"); + + Teuchos::ParameterList pl; + pl.set("Basis", basis_qedge1); + pl.set("DOF Names", tangent_names); + pl.set("Indexer Names", names); + + pl.set("Global Data Key", "Tangent Container"); + + 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); + ped.gedc->addDataObject("Tangent Container", tangentContainers[0]); + std::vector params; + params.push_back("param1"); + params.push_back("param2"); + params.push_back("param3"); + Teuchos::RCP activeParams = + Teuchos::rcp(new panzer::ParameterList_GlobalEvaluationData(params)); + ped.gedc->addDataObject("PARAMETER_NAMES",activeParams); + auto paramContainer1 = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F)); + ped.gedc->addDataObject("param1",paramContainer1->getGhostedLOC()); + auto paramContainer2 = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F)); + ped.gedc->addDataObject("param2",paramContainer2->getGhostedLOC()); + auto paramContainer3 = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F)); + ped.gedc->addDataObject("param3",paramContainer3->getGhostedLOC()); + + 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); + + // test Tangent fields + { + Teuchos::ArrayRCP data; + Teuchos::RCP> f_vec = t_loc->get_f_th(); + + // 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)->getLocalData(Teuchos::ptrFromRef(data)); + for (size_type i = 0; i < data.size(); i++) + { + double target = 123.0 + myRank; + TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); + } + } + { + // now check the tangent values + Teuchos::ArrayRCP data1, data2; + Teuchos::RCP> f1_vec = Teuchos::rcp_dynamic_cast(paramContainer1->getGhostedLOC())->get_f_th(); + Teuchos::rcp_dynamic_cast>(f1_vec)->getLocalData(Teuchos::ptrFromRef(data1)); + + Teuchos::RCP> f2_vec = Teuchos::rcp_dynamic_cast(paramContainer2->getGhostedLOC())->get_f_th(); + Teuchos::rcp_dynamic_cast>(f2_vec)->getLocalData(Teuchos::ptrFromRef(data2)); + + for (size_type i = 0; i < data1.size(); ++i) + { + std::cout << "DATA :: " << data1[i] << " " << data2[i] << std::endl; + } + } + + PHX::MDField + fieldData1_q1(fieldName1_q1, basis_q1->functional); + PHX::MDField + fieldData2_q1(fieldName2_q1, basis_q1->functional); + + fm.getFieldData(fieldData1_q1); + fm.getFieldData(fieldData2_q1); + + auto fieldData1_q1_h = Kokkos::create_mirror_view(fieldData1_q1.get_static_view()); + auto fieldData2_q1_h = Kokkos::create_mirror_view(fieldData2_q1.get_static_view()); + Kokkos::deep_copy(fieldData1_q1_h, fieldData1_q1.get_static_view()); + Kokkos::deep_copy(fieldData2_q1_h, fieldData2_q1.get_static_view()); + for (unsigned int cell = 0; cell < fieldData1_q1.extent(0); ++cell) + { + for (unsigned int pt = 0; pt < fieldData1_q1.extent(1); pt++) + { + std::cout << fieldData1_q1_h(cell, pt).dx(0) << " " << fieldData2_q1_h(cell,pt).dx(0) << std::endl; + } + } + } + Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName) { Teuchos::RCP topo = diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 7639dc582391..4c8c57d914a2 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -206,11 +206,12 @@ preEvaluate(typename TRAITS::PreEvalData d) dfdpFieldsVoV_.initialize("ScatterResidual_Tpetra::dfdpFieldsVoV_",activeParameters.size()); for(std::size_t i=0;if dfdp? RCP vec = rcp_dynamic_cast(d.gedc->getDataObject(activeParameters[i]),true)->get_f(); auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access::ReadWrite); - // TODO BWR not sure this will be checked by tests... + // TODO BWR not sure this will be checked by tests... Always true?? TEUCHOS_ASSERT(dfdp_view.extent_int(1)==1); dfdpFieldsVoV_.addView(Kokkos::subview(dfdp_view,Kokkos::ALL(),0),i); } @@ -409,15 +410,14 @@ class ScatterResidual_Tangent_Functor { Kokkos::View lids; // local indices for unknowns. PHX::View offsets; // how to get a particular field FieldType field; + double num_derivs; - PHX::ViewOfViews<1,PHX::View> dfdp_fields; // tangent fields + Kokkos::View*> dfdp_fields; // tangent fields KOKKOS_INLINE_FUNCTION void operator()(const unsigned int cell) const { - const std::size_t numDerivs = PHX::getFadSize(field.get_static_view()); - const auto tangents = dfdp_fields.getViewDevice(); // loop over the basis functions (currently they are nodes) for(std::size_t basis=0; basis < offsets.extent(0); basis++) { typename FieldType::array_type::reference_type scatterField = field(cell,basis); @@ -429,8 +429,9 @@ class ScatterResidual_Tangent_Functor { Kokkos::atomic_add(&r_data(lid,0), scatterField.val()); // loop over the tangents - for(int i_param=0; i_param Date: Fri, 8 Nov 2024 12:25:13 -0700 Subject: [PATCH 04/11] Panzer Tangents :: Fix unit test Signed-off-by: Bryan Reuter --- .../tpetra_scatter_residual.cpp | 136 ++++++++++-------- .../Panzer_GatherSolution_Tpetra_impl.hpp | 4 +- .../Panzer_ScatterResidual_Tpetra_impl.hpp | 7 +- 3 files changed, 80 insertions(+), 67 deletions(-) diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp index 3c63810a0cd5..a360232cb795 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp @@ -65,8 +65,6 @@ namespace panzer void testInitialization(const Teuchos::RCP &ipb); Teuchos::RCP buildMesh(int elemX, int elemY); - // TODO BWR NO TANGENT COVERAGE HERE OR IN OTHER TEST - TEUCHOS_UNIT_TEST(assembly, scatter_solution_residual) { @@ -454,7 +452,6 @@ namespace panzer TEUCHOS_UNIT_TEST(assembly, scatter_solution_tangent) { - // TODO BWR not checking the tangent calculation, just the residual part!!! #ifdef HAVE_MPI Teuchos::RCP> tComm = Teuchos::rcp(new Teuchos::MpiComm(MPI_COMM_WORLD)); @@ -468,6 +465,7 @@ namespace panzer const std::string fieldName1_q1 = "U"; const std::string fieldName2_q1 = "V"; const std::string fieldName_qedge1 = "B"; + const std::size_t numParams = 3; Teuchos::RCP mesh = buildMesh(2, 2); @@ -523,17 +521,20 @@ namespace panzer using LOCPair = panzer::LOCPair_GlobalEvaluationData; using Teuchos::rcp_dynamic_cast; - auto locPair = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X)); + // generate tangent data + for (std::size_t i=0;i(locPair->getGlobalLOC()); - Teuchos::RCP> global_x_vec = global_t_loc->get_x_th(); - Thyra::assign(global_x_vec.ptr(), 0.123 + myRank); + auto global_t_loc = rcp_dynamic_cast(locPair->getGlobalLOC()); + Teuchos::RCP> global_x_vec = global_t_loc->get_x_th(); + Thyra::assign(global_x_vec.ptr(), 0.123 + myRank + i); - auto ghosted_t_loc = rcp_dynamic_cast(locPair->getGhostedLOC()); - Teuchos::RCP> ghosted_x_vec = ghosted_t_loc->get_x_th(); - Thyra::assign(ghosted_x_vec.ptr(), 0.123 + myRank); + auto ghosted_t_loc = rcp_dynamic_cast(locPair->getGhostedLOC()); + Teuchos::RCP> ghosted_x_vec = ghosted_t_loc->get_x_th(); + Thyra::assign(ghosted_x_vec.ptr(), 0.123 + myRank + i); - tangentContainers.push_back(locPair); + tangentContainers.push_back(locPair); + } // setup field manager, add evaluator under test ///////////////////////////////////////////////////////////// @@ -541,7 +542,7 @@ namespace panzer PHX::FieldManager fm; std::vector derivative_dimensions; - derivative_dimensions.push_back(1); + derivative_dimensions.push_back(numParams); fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); std::string resName = ""; @@ -551,6 +552,13 @@ namespace panzer names_map->insert(std::make_pair(fieldName2_q1, resName + fieldName2_q1)); names_map->insert(std::make_pair(fieldName_qedge1, resName + fieldName_qedge1)); + // Guide to what's happening in this test + // 1) U,V,B gathers request tangent fields so those are first gathered with gatherTangent + // 2) When U,V,B gathers occur, the tangents are stored as derivatives. The first tangent is the first derivative and so on. + // 3) U,V,B are overloaded to also be the residual fields we are scattering + // 4) This means their derivatives, e.g., U.dx(i) are considered to be dfdp_i -> derivatives of the residuals w.r.t. parameters. + // 5) dfdp_i = f.dx(i) so the number of tangents is the number of params + // evaluators under test { using Teuchos::RCP; @@ -606,30 +614,43 @@ namespace panzer pl.set("Indexer Names", names); Teuchos::RCP>> tangent_names = Teuchos::rcp(new std::vector>(2)); - (*tangent_names)[0].push_back(fieldName1_q1 + " Tangent"); - (*tangent_names)[1].push_back(fieldName2_q1 + " Tangent"); + for (int 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); - tangent_names->push_back(fieldName1_q1 + " Tangent"); - tangent_names->push_back(fieldName2_q1 + " Tangent"); + { + 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); - pl.set("Global Data Key", "Tangent Container"); + std::stringstream ss; + ss << "Tangent Container " << i; + pl.set("Global Data Key", ss.str()); Teuchos::RCP> evaluator = lof->buildGatherTangent(pl); @@ -648,34 +669,44 @@ namespace panzer pl.set("Indexer Names", names); Teuchos::RCP>> tangent_names = Teuchos::rcp(new std::vector>(1)); - (*tangent_names)[0].push_back(fieldName_qedge1 + " Tangent"); + for (int 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); - tangent_names->push_back(fieldName_qedge1 + " Tangent"); + { + 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); - pl.set("Global Data Key", "Tangent Container"); + 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; @@ -684,20 +715,24 @@ namespace panzer panzer::Traits::PED ped; ped.gedc->addDataObject("Solution Gather Container", loc); ped.gedc->addDataObject("Residual Scatter Container", loc); - ped.gedc->addDataObject("Tangent Container", tangentContainers[0]); + for (size_t i=0; iaddDataObject(ss.str(), tangentContainers[i]); + } std::vector params; - params.push_back("param1"); - params.push_back("param2"); - params.push_back("param3"); + 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); - auto paramContainer1 = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F)); - ped.gedc->addDataObject("param1",paramContainer1->getGhostedLOC()); - auto paramContainer2 = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F)); - ped.gedc->addDataObject("param2",paramContainer2->getGhostedLOC()); - auto paramContainer3 = Teuchos::rcp(new LOCPair(t_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F)); - ped.gedc->addDataObject("param3",paramContainer3->getGhostedLOC()); fm.preEvaluate(ped); @@ -728,40 +763,19 @@ namespace panzer TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); } } + for (std::size_t i=0; i data1, data2; - Teuchos::RCP> f1_vec = Teuchos::rcp_dynamic_cast(paramContainer1->getGhostedLOC())->get_f_th(); - Teuchos::rcp_dynamic_cast>(f1_vec)->getLocalData(Teuchos::ptrFromRef(data1)); - - Teuchos::RCP> f2_vec = Teuchos::rcp_dynamic_cast(paramContainer2->getGhostedLOC())->get_f_th(); - Teuchos::rcp_dynamic_cast>(f2_vec)->getLocalData(Teuchos::ptrFromRef(data2)); + Teuchos::ArrayRCP data; + Teuchos::RCP> f_vec = Teuchos::rcp_dynamic_cast(paramContainers[i]->getGhostedLOC())->get_f_th(); + Teuchos::rcp_dynamic_cast>(f_vec)->getLocalData(Teuchos::ptrFromRef(data)); - for (size_type i = 0; i < data1.size(); ++i) + for (std::size_t j = 0; j < data.size(); ++j) { - std::cout << "DATA :: " << data1[i] << " " << data2[i] << std::endl; + const double target = .123 + myRank + i; + TEST_ASSERT(data[j] == target || data[j] == 2.0 * target); } } - - PHX::MDField - fieldData1_q1(fieldName1_q1, basis_q1->functional); - PHX::MDField - fieldData2_q1(fieldName2_q1, basis_q1->functional); - - fm.getFieldData(fieldData1_q1); - fm.getFieldData(fieldData2_q1); - - auto fieldData1_q1_h = Kokkos::create_mirror_view(fieldData1_q1.get_static_view()); - auto fieldData2_q1_h = Kokkos::create_mirror_view(fieldData2_q1.get_static_view()); - Kokkos::deep_copy(fieldData1_q1_h, fieldData1_q1.get_static_view()); - Kokkos::deep_copy(fieldData2_q1_h, fieldData2_q1.get_static_view()); - for (unsigned int cell = 0; cell < fieldData1_q1.extent(0); ++cell) - { - for (unsigned int pt = 0; pt < fieldData1_q1.extent(1); pt++) - { - std::cout << fieldData1_q1_h(cell, pt).dx(0) << " " << fieldData2_q1_h(cell,pt).dx(0) << std::endl; - } - } } Teuchos::RCP buildBasis(std::size_t worksetSize, const std::string &basisName) diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_Tpetra_impl.hpp index bdff0f73574e..21be837adf3e 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_Tpetra_impl.hpp @@ -263,8 +263,8 @@ postRegistrationSetup(typename TRAITS::SetupData /* d */, } Kokkos::deep_copy(tangentInnerVectorSizes_,tangentInnerVectorSizes_host); - gatherFieldsVoV_.initialize("GatherSolution_Teptra::gatherFieldsVoV_",gatherFields_.size()); - tangentFieldsVoV_.initialize("GatherSolution_Teptra::tangentFieldsVoV_",gatherFields_.size(),inner_vector_max_size); + gatherFieldsVoV_.initialize("GatherSolution_Tpetra::gatherFieldsVoV_",gatherFields_.size()); + tangentFieldsVoV_.initialize("GatherSolution_Tpetra::tangentFieldsVoV_",gatherFields_.size(),inner_vector_max_size); for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) { const std::string& fieldName = indexerNames_[fd]; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 4c8c57d914a2..94be35b0d188 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -410,7 +410,7 @@ class ScatterResidual_Tangent_Functor { Kokkos::View lids; // local indices for unknowns. PHX::View offsets; // how to get a particular field FieldType field; - double num_derivs; + double num_params; Kokkos::View*> dfdp_fields; // tangent fields @@ -429,8 +429,7 @@ class ScatterResidual_Tangent_Functor { Kokkos::atomic_add(&r_data(lid,0), scatterField.val()); // loop over the tangents - // TODO BWR HERE WE ARE HITTING WRONG SLOT WITH TEST!! - for(int i_param=0; i_param Date: Fri, 8 Nov 2024 13:50:52 -0700 Subject: [PATCH 05/11] Panzer Tangents :: wsign error Signed-off-by: Bryan Reuter --- .../test/evaluator_tests/tpetra_scatter_residual.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp index a360232cb795..0ae80e3d9b50 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp @@ -763,14 +763,14 @@ namespace panzer TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); } } - for (std::size_t i=0; i data; Teuchos::RCP> f_vec = Teuchos::rcp_dynamic_cast(paramContainers[i]->getGhostedLOC())->get_f_th(); Teuchos::rcp_dynamic_cast>(f_vec)->getLocalData(Teuchos::ptrFromRef(data)); - for (std::size_t j = 0; j < data.size(); ++j) + for (size_type j = 0; j < data.size(); ++j) { const double target = .123 + myRank + i; TEST_ASSERT(data[j] == target || data[j] == 2.0 * target); From a6be01e72247c0915a5621bbdad026104fbb7e99 Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Fri, 8 Nov 2024 14:22:22 -0700 Subject: [PATCH 06/11] Panzer Tangents :: more wsign error Signed-off-by: Bryan Reuter --- .../test/evaluator_tests/tpetra_scatter_residual.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp index 0ae80e3d9b50..9fc68bab53ac 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp @@ -614,7 +614,7 @@ namespace panzer pl.set("Indexer Names", names); Teuchos::RCP>> tangent_names = Teuchos::rcp(new std::vector>(2)); - for (int i = 0; i < numParams; ++i) + for (std::size_t i = 0; i < numParams; ++i) { std::stringstream ss1, ss2; ss1 << fieldName1_q1 << " Tangent " << i; @@ -669,7 +669,7 @@ namespace panzer pl.set("Indexer Names", names); Teuchos::RCP>> tangent_names = Teuchos::rcp(new std::vector>(1)); - for (int i = 0; i < numParams; ++i) + for (std::size_t i = 0; i < numParams; ++i) { std::stringstream ss; ss << fieldName_qedge1 << " Tangent " << i; @@ -763,7 +763,7 @@ namespace panzer TEST_ASSERT(data[i] == target || data[i] == 2.0 * target); } } - for (int i=0; i data; From 53dff92077d478445e2cccfef5bc77e567cff8f9 Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Fri, 8 Nov 2024 16:02:47 -0700 Subject: [PATCH 07/11] Panzer Tangents :: Hopefully fix cuda compile error Signed-off-by: Bryan Reuter --- .../src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 94be35b0d188..ec259b4f3e65 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -412,7 +412,7 @@ class ScatterResidual_Tangent_Functor { FieldType field; double num_params; - Kokkos::View*> dfdp_fields; // tangent fields + Kokkos::View*> dfdp_fields; // tangent fields KOKKOS_INLINE_FUNCTION void operator()(const unsigned int cell) const From 2d788b000abc44f53e06ada34ff6e921e64f8127 Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Sat, 9 Nov 2024 16:06:25 -0700 Subject: [PATCH 08/11] Panzer Tangents :: Compiling on cuda Signed-off-by: Bryan Reuter --- .../src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index ec259b4f3e65..8c9afbb8fbd3 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -407,7 +407,7 @@ class ScatterResidual_Tangent_Functor { bool fillResidual; Kokkos::View r_data; - Kokkos::View lids; // local indices for unknowns. + Kokkos::View lids; // local indices for unknowns. PHX::View offsets; // how to get a particular field FieldType field; double num_params; From bdd852489965bcbdcf927e7374ebf69b40304f2d Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Mon, 11 Nov 2024 15:44:42 -0700 Subject: [PATCH 09/11] Panzer Tangents :: Release FM in unit test and address Roger's comments Signed-off-by: Bryan Reuter --- .../tpetra_scatter_residual.cpp | 30 ++++++++++--------- .../Panzer_ScatterResidual_Tpetra_decl.hpp | 2 +- .../Panzer_ScatterResidual_Tpetra_impl.hpp | 15 ++++------ 3 files changed, 22 insertions(+), 25 deletions(-) diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp index 9fc68bab53ac..1e86353a14bc 100644 --- a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp +++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_scatter_residual.cpp @@ -539,11 +539,11 @@ namespace panzer // setup field manager, add evaluator under test ///////////////////////////////////////////////////////////// - PHX::FieldManager fm; + auto fm = Teuchos::rcp(new PHX::FieldManager); std::vector derivative_dimensions; derivative_dimensions.push_back(numParams); - fm.setKokkosExtendedDataTypeDimensions(derivative_dimensions); + fm->setKokkosExtendedDataTypeDimensions(derivative_dimensions); std::string resName = ""; Teuchos::RCP> names_map = @@ -577,8 +577,8 @@ namespace panzer TEST_EQUALITY(evaluator->evaluatedFields().size(), 1); - fm.registerEvaluator(evaluator); - fm.requireField(*evaluator->evaluatedFields()[0]); + fm->registerEvaluator(evaluator); + fm->requireField(*evaluator->evaluatedFields()[0]); } { using Teuchos::RCP; @@ -596,8 +596,8 @@ namespace panzer TEST_EQUALITY(evaluator->evaluatedFields().size(), 1); - fm.registerEvaluator(evaluator); - fm.requireField(*evaluator->evaluatedFields()[0]); + fm->registerEvaluator(evaluator); + fm->requireField(*evaluator->evaluatedFields()[0]); } // support evaluators @@ -626,7 +626,7 @@ namespace panzer Teuchos::RCP> evaluator = lof->buildGather(pl); - fm.registerEvaluator(evaluator); + fm->registerEvaluator(evaluator); } for (std::size_t i = 0; i < numParams; ++i) { using Teuchos::RCP; @@ -655,7 +655,7 @@ namespace panzer Teuchos::RCP> evaluator = lof->buildGatherTangent(pl); - fm.registerEvaluator(evaluator); + fm->registerEvaluator(evaluator); } { using Teuchos::RCP; @@ -679,7 +679,7 @@ namespace panzer Teuchos::RCP> evaluator = lof->buildGather(pl); - fm.registerEvaluator(evaluator); + fm->registerEvaluator(evaluator); } for (std::size_t i = 0; i < numParams; ++i) { using Teuchos::RCP; @@ -705,12 +705,12 @@ namespace panzer Teuchos::RCP> evaluator = lof->buildGatherTangent(pl); - fm.registerEvaluator(evaluator); + fm->registerEvaluator(evaluator); } panzer::Traits::SD sd; sd.worksets_ = work_sets; - fm.postRegistrationSetup(sd); + fm->postRegistrationSetup(sd); panzer::Traits::PED ped; ped.gedc->addDataObject("Solution Gather Container", loc); @@ -734,7 +734,7 @@ namespace panzer Teuchos::rcp(new panzer::ParameterList_GlobalEvaluationData(params)); ped.gedc->addDataObject("PARAMETER_NAMES",activeParams); - fm.preEvaluate(ped); + fm->preEvaluate(ped); // run tests ///////////////////////////////////////////////////////////// @@ -745,8 +745,10 @@ namespace panzer workset.time = 0.0; workset.evaluate_transient_terms = false; - fm.evaluateFields(workset); - fm.postEvaluate(0); + fm->evaluateFields(workset); + fm->postEvaluate(0); + + fm = Teuchos::null; // test Tangent fields { diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp index 04609a23a076..ad06fbe8f2d7 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_decl.hpp @@ -158,7 +158,7 @@ class ScatterResidual_Tpetra Teuchos::RCP > tpetraContainer_; /// Storage for the tangent data - PHX::ViewOfViews<1,PHX::View> dfdpFieldsVoV_; + PHX::ViewOfViews<1,Kokkos::View> dfdpFieldsVoV_; PHX::View scratch_lids_; std::vector > scratch_offsets_; diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 8c9afbb8fbd3..09420fcfd4da 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -24,7 +24,6 @@ #include "Panzer_GlobalEvaluationDataContainer.hpp" #include "Phalanx_DataLayout_MDALayout.hpp" -#include "Phalanx_Scratch_Utilities.hpp" #include "Teuchos_FancyOStream.hpp" @@ -211,14 +210,11 @@ preEvaluate(typename TRAITS::PreEvalData d) rcp_dynamic_cast(d.gedc->getDataObject(activeParameters[i]),true)->get_f(); auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access::ReadWrite); - // TODO BWR not sure this will be checked by tests... Always true?? - TEUCHOS_ASSERT(dfdp_view.extent_int(1)==1); - dfdpFieldsVoV_.addView(Kokkos::subview(dfdp_view,Kokkos::ALL(),0),i); + dfdpFieldsVoV_.addView(dfdp_view,i); } dfdpFieldsVoV_.syncHostToDevice(); - // TODO BWR DO WE WANT TO BE ABLE TO FILL RESIDUAL TOO? // extract linear object container tpetraContainer_ = Teuchos::rcp_dynamic_cast(d.gedc->getDataObject(globalDataKey_)); @@ -412,7 +408,7 @@ class ScatterResidual_Tangent_Functor { FieldType field; double num_params; - Kokkos::View*> dfdp_fields; // tangent fields + Kokkos::View*> dfdp_fields; // tangent fields KOKKOS_INLINE_FUNCTION void operator()(const unsigned int cell) const @@ -430,7 +426,7 @@ class ScatterResidual_Tangent_Functor { // loop over the tangents for(int i_param=0; i_param functor; functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite); functor.lids = scratch_lids_; + functor.dfdp_fields = dfdpFieldsVoV_.getViewDevice(); // for each field, do a parallel for loop for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) { functor.offsets = scratch_offsets_[fieldIndex]; functor.field = scatterFields_[fieldIndex]; - functor.dfdp_fields = dfdpFieldsVoV_.getViewDevice(); - // TODO BWR is there a better way? Is this correct? - functor.num_params = PHX::getFadSize(scatterFields_[fieldIndex].get_static_view())-1; + functor.num_params = Kokkos::dimension_scalar(scatterFields_[fieldIndex].get_view())-1; Kokkos::parallel_for(workset.num_cells,functor); } From 8ff1aeb1a5c31fa42b5b62bf82ce33f7e4f7b184 Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Tue, 12 Nov 2024 14:38:28 -0700 Subject: [PATCH 10/11] Panzer Tangents :: Final clean up Signed-off-by: Bryan Reuter --- ...zer_ScatterResidual_BlockedTpetra_impl.hpp | 82 +++++++++++++++++++ .../Panzer_ScatterResidual_Tpetra_impl.hpp | 1 - 2 files changed, 82 insertions(+), 1 deletion(-) 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 3e8633c7f3ef..8cf1e9529d5a 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 @@ -472,6 +472,88 @@ 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); +} + // ********************************************************************** #endif diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp index 09420fcfd4da..c3655218f381 100644 --- a/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp +++ b/packages/panzer/disc-fe/src/evaluators/Panzer_ScatterResidual_Tpetra_impl.hpp @@ -205,7 +205,6 @@ preEvaluate(typename TRAITS::PreEvalData d) dfdpFieldsVoV_.initialize("ScatterResidual_Tpetra::dfdpFieldsVoV_",activeParameters.size()); for(std::size_t i=0;if dfdp? RCP vec = rcp_dynamic_cast(d.gedc->getDataObject(activeParameters[i]),true)->get_f(); auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access::ReadWrite); From 7c544a6821704d0296967414ea25c51613aa319b Mon Sep 17 00:00:00 2001 From: Bryan Reuter Date: Tue, 12 Nov 2024 15:03:24 -0700 Subject: [PATCH 11/11] Panzer Tangents :: Remove blocked prototype Signed-off-by: Bryan Reuter --- ...zer_ScatterResidual_BlockedTpetra_impl.hpp | 84 ------------------- 1 file changed, 84 deletions(-) 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 8cf1e9529d5a..743d8b30cef7 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 @@ -472,88 +472,4 @@ 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); -} - -// ********************************************************************** - #endif