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