Skip to content

Commit

Permalink
Panzer tangent scatter blocked tpetra (trilinos#13601)
Browse files Browse the repository at this point in the history 
* Panzer Tangents :: Scatter blocked Tpetra updated with unit test

---------

Signed-off-by: Bryan Reuter <[email protected]>
  • Loading branch information
@reuterb
reuterb committed Nov 26, 2024
1 parent cb93f35 commit 49a91f1
Show file tree
Hide file tree
Showing 5 changed files with 716 additions and 75 deletions.
363 changes: 363 additions & 0 deletions packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_scatter_residual.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,9 @@ using Teuchos::rcp;
#include "Panzer_GatherOrientation.hpp"
#include "Panzer_ScatterResidual_BlockedTpetra.hpp"
#include "Panzer_GatherSolution_BlockedTpetra.hpp"
#include "Panzer_LOCPair_GlobalEvaluationData.hpp"
#include "Panzer_GlobalEvaluationDataContainer.hpp"
#include "Panzer_ParameterList_GlobalEvaluationData.hpp"

#include "Panzer_STK_Version.hpp"
#include "PanzerAdaptersSTK_config.hpp"
Expand Down Expand Up @@ -494,6 +496,367 @@ namespace panzer
}
}

TEUCHOS_UNIT_TEST(block_assembly, scatter_solution_tangent)
{

#ifdef HAVE_MPI
Teuchos::RCP<Teuchos::MpiComm<int>> tComm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));
#else
NOPE_PANZER_DOESNT_SUPPORT_SERIAL
#endif

int myRank = tComm->getRank();

const std::size_t workset_size = 4;
const std::string fieldName1_q1 = "U";
const std::string fieldName2_q1 = "V";
const std::string fieldName_qedge1 = "B";
const std::size_t numParams = 2;

Teuchos::RCP<panzer_stk::STK_Interface> mesh = buildMesh(2, 2);

// build input physics block
Teuchos::RCP<panzer::PureBasis> basis_q1 = buildBasis(workset_size, "Q1");
Teuchos::RCP<panzer::PureBasis> basis_qedge1 = buildBasis(workset_size, "QEdge1");

Teuchos::RCP<Teuchos::ParameterList> ipb = Teuchos::parameterList();
testInitialization(ipb);

const int default_int_order = 1;
std::string eBlockID = "eblock-0_0";
Teuchos::RCP<user_app::MyFactory> eqset_factory = Teuchos::rcp(new user_app::MyFactory);
panzer::CellData cellData(workset_size, mesh->getCellTopology("eblock-0_0"));
Teuchos::RCP<panzer::GlobalData> gd = panzer::createGlobalData();
Teuchos::RCP<panzer::PhysicsBlock> physicsBlock =
Teuchos::rcp(new PhysicsBlock(ipb, eBlockID, default_int_order, cellData, eqset_factory, gd, false));

Teuchos::RCP<std::vector<panzer::Workset>> 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<panzer::ConnManager> conn_manager = Teuchos::rcp(new panzer_stk::STKConnManager(mesh));
RCP<panzer::BlockedDOFManager> 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<std::vector<std::string>> 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<TpetraBlockedLinObjFactoryType> bt_lof = Teuchos::rcp(new TpetraBlockedLinObjFactoryType(tComm.getConst(), dofManager));
Teuchos::RCP<LinearObjFactory<panzer::Traits>> lof = bt_lof;
Teuchos::RCP<LinearObjContainer> loc = bt_lof->buildGhostedLinearObjContainer();
bt_lof->initializeGhostedContainer(LinearObjContainer::X | LinearObjContainer::F, *loc);
loc->initialize();

Teuchos::RCP<TpetraBlockedLinObjContainerType> b_loc = Teuchos::rcp_dynamic_cast<TpetraBlockedLinObjContainerType>(loc);

Teuchos::RCP<Thyra::ProductVectorBase<double>> p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double>>(b_loc->get_x());
Thyra::assign(p_vec->getNonconstVectorBlock(0).ptr(), 123.0 + myRank);
Thyra::assign(p_vec->getNonconstVectorBlock(1).ptr(), 456.0 + myRank);
Thyra::assign(p_vec->getNonconstVectorBlock(2).ptr(), 789.0 + myRank);

std::vector<Teuchos::RCP<GlobalEvaluationData>> tangentContainers;

using LOCPair = panzer::LOCPair_GlobalEvaluationData;
using Teuchos::rcp_dynamic_cast;

// generate tangent data
for (std::size_t i=0;i<numParams; ++i){
auto locPair = Teuchos::rcp(new LOCPair(bt_lof, panzer::LinearObjContainer::X));

auto global_bt_loc = rcp_dynamic_cast<TpetraBlockedLinObjContainerType>(locPair->getGlobalLOC());
Teuchos::RCP<Thyra::ProductVectorBase<double>> global_p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double>>(global_bt_loc->get_x());
Thyra::assign(global_p_vec->getNonconstVectorBlock(0).ptr(), 0.123 + myRank + i);
Thyra::assign(global_p_vec->getNonconstVectorBlock(1).ptr(), 0.456 + myRank + i);
Thyra::assign(global_p_vec->getNonconstVectorBlock(2).ptr(), 0.789 + myRank + i);

auto ghosted_bt_loc = rcp_dynamic_cast<TpetraBlockedLinObjContainerType>(locPair->getGhostedLOC());
Teuchos::RCP<Thyra::ProductVectorBase<double>> ghosted_p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double>>(ghosted_bt_loc->get_x());
Thyra::assign(ghosted_p_vec->getNonconstVectorBlock(0).ptr(), 0.123 + myRank + i);
Thyra::assign(ghosted_p_vec->getNonconstVectorBlock(1).ptr(), 0.456 + myRank + i);
Thyra::assign(ghosted_p_vec->getNonconstVectorBlock(2).ptr(), 0.789 + myRank + i);

tangentContainers.push_back(locPair);
}

// setup field manager, add evaluator under test
/////////////////////////////////////////////////////////////

auto fm = Teuchos::rcp(new PHX::FieldManager<panzer::Traits>);

std::vector<PHX::index_size_type> derivative_dimensions;
derivative_dimensions.push_back(numParams);
fm->setKokkosExtendedDataTypeDimensions<panzer::Traits::Tangent>(derivative_dimensions);

std::string resName = "";
Teuchos::RCP<std::map<std::string, std::string>> names_map =
Teuchos::rcp(new std::map<std::string, std::string>);
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<std::vector<std::string>> names = rcp(new std::vector<std::string>);
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<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildScatter<panzer::Traits::Tangent>(pl);

TEST_EQUALITY(evaluator->evaluatedFields().size(), 1);

fm->registerEvaluator<panzer::Traits::Tangent>(evaluator);
fm->requireField<panzer::Traits::Tangent>(*evaluator->evaluatedFields()[0]);
}
{
using Teuchos::RCP;
using Teuchos::rcp;
RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
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<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildScatter<panzer::Traits::Tangent>(pl);

TEST_EQUALITY(evaluator->evaluatedFields().size(), 1);

fm->registerEvaluator<panzer::Traits::Tangent>(evaluator);
fm->requireField<panzer::Traits::Tangent>(*evaluator->evaluatedFields()[0]);
}

// support evaluators
{
using Teuchos::RCP;
using Teuchos::rcp;
RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
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<std::vector<std::vector<std::string>>> tangent_names =
Teuchos::rcp(new std::vector<std::vector<std::string>>(2));
for (std::size_t i = 0; i < numParams; ++i)
{
std::stringstream ss1, ss2;
ss1 << fieldName1_q1 << " Tangent " << i;
ss2 << fieldName2_q1 << " Tangent " << i;
(*tangent_names)[0].push_back(ss1.str());
(*tangent_names)[1].push_back(ss2.str());
}
pl.set("Tangent Names", tangent_names);

Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Tangent>(pl);

fm->registerEvaluator<panzer::Traits::Tangent>(evaluator);
}
for (std::size_t i = 0; i < numParams; ++i) {
using Teuchos::RCP;
using Teuchos::rcp;
RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
RCP<std::vector<std::string>> tangent_names = rcp(new std::vector<std::string>);
names->push_back(fieldName1_q1);
names->push_back(fieldName2_q1);
{
std::stringstream ss1, ss2;
ss1 << fieldName1_q1 << " Tangent " << i;
ss2 << fieldName2_q1 << " Tangent " << i;
tangent_names->push_back(ss1.str());
tangent_names->push_back(ss2.str());
}

Teuchos::ParameterList pl;
pl.set("Basis", basis_q1);
pl.set("DOF Names", tangent_names);
pl.set("Indexer Names", names);

std::stringstream ss;
ss << "Tangent Container " << i;
pl.set("Global Data Key", ss.str());

Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator =
lof->buildGatherTangent<panzer::Traits::Tangent>(pl);

fm->registerEvaluator<panzer::Traits::Tangent>(evaluator);
}
{
using Teuchos::RCP;
using Teuchos::rcp;
RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
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<std::vector<std::vector<std::string>>> tangent_names =
Teuchos::rcp(new std::vector<std::vector<std::string>>(1));
for (std::size_t i = 0; i < numParams; ++i)
{
std::stringstream ss;
ss << fieldName_qedge1 << " Tangent " << i;
(*tangent_names)[0].push_back(ss.str());
}
pl.set("Tangent Names", tangent_names);

Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Tangent>(pl);

fm->registerEvaluator<panzer::Traits::Tangent>(evaluator);
}
for (std::size_t i = 0; i < numParams; ++i) {
using Teuchos::RCP;
using Teuchos::rcp;
RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
RCP<std::vector<std::string>> tangent_names = rcp(new std::vector<std::string>);
names->push_back(fieldName_qedge1);
{
std::stringstream ss;
ss << fieldName_qedge1 << " Tangent " << i;
tangent_names->push_back(ss.str());
}

Teuchos::ParameterList pl;
pl.set("Basis", basis_qedge1);
pl.set("DOF Names", tangent_names);
pl.set("Indexer Names", names);

std::stringstream ss;
ss << "Tangent Container " << i;
pl.set("Global Data Key", ss.str());

Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator =
lof->buildGatherTangent<panzer::Traits::Tangent>(pl);

fm->registerEvaluator<panzer::Traits::Tangent>(evaluator);
}

panzer::Traits::SD sd;
sd.worksets_ = work_sets;

fm->postRegistrationSetup(sd);

panzer::Traits::PED ped;
ped.gedc->addDataObject("Solution Gather Container", loc);
ped.gedc->addDataObject("Residual Scatter Container", loc);
for (size_t i=0; i<numParams; ++i) {
std::stringstream ss;
ss << "Tangent Container " << i;
ped.gedc->addDataObject(ss.str(), tangentContainers[i]);
}
std::vector<std::string> params;
std::vector<Teuchos::RCP<LOCPair>> paramContainers;
for (std::size_t i = 0; i<numParams; ++i) {
std::stringstream ss;
ss << "param" << i;
params.push_back(ss.str());
auto paramContainer = Teuchos::rcp(new LOCPair(bt_lof, panzer::LinearObjContainer::X | panzer::LinearObjContainer::F));
ped.gedc->addDataObject(ss.str(),paramContainer->getGhostedLOC());
paramContainers.push_back(paramContainer);
}
Teuchos::RCP<panzer::GlobalEvaluationData> activeParams =
Teuchos::rcp(new panzer::ParameterList_GlobalEvaluationData(params));
ped.gedc->addDataObject("PARAMETER_NAMES",activeParams);
fm->preEvaluate<panzer::Traits::Tangent>(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<panzer::Traits::Tangent>(workset);
fm->postEvaluate<panzer::Traits::Tangent>(0);

fm = Teuchos::null;

// test Tangent fields
Teuchos::ArrayRCP<const double> data;
Teuchos::RCP<const Thyra::ProductVectorBase<double>> f_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double>>(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<const Thyra::SpmdVectorBase<double>>(f_vec->getVectorBlock(0))->getLocalData(Teuchos::ptrFromRef(data));
TEST_EQUALITY(static_cast<size_t>(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<const Thyra::SpmdVectorBase<double>>(f_vec->getVectorBlock(1))->getLocalData(Teuchos::ptrFromRef(data));
TEST_EQUALITY(static_cast<size_t>(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<const Thyra::SpmdVectorBase<double>>(f_vec->getVectorBlock(2))->getLocalData(Teuchos::ptrFromRef(data));
TEST_EQUALITY(static_cast<size_t>(data.size()), b_loc->getMapForBlock(2)->getLocalNumElements());
for (size_type i = 0; i < data.size(); i++)
{
double target = 789.0 + myRank;
TEST_ASSERT(data[i] == target || data[i] == 2.0 * target);
}

for (std::size_t i=0; i<numParams; ++i) {
Teuchos::RCP<const Thyra::ProductVectorBase<double>> param_f_vec =
Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double>>(
Teuchos::rcp_dynamic_cast<TpetraBlockedLinObjContainerType>(paramContainers[i]->getGhostedLOC())->get_f());

Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorBase<double>>(param_f_vec->getVectorBlock(0))->getLocalData(Teuchos::ptrFromRef(data));
TEST_EQUALITY(static_cast<size_t>(data.size()), b_loc->getMapForBlock(0)->getLocalNumElements());
for (size_type j = 0; j < data.size(); j++)
{
double target = .123 + myRank + i;
TEST_ASSERT(data[j] == target || data[j] == 2.0 * target);
}
Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorBase<double>>(param_f_vec->getVectorBlock(1))->getLocalData(Teuchos::ptrFromRef(data));
TEST_EQUALITY(static_cast<size_t>(data.size()), b_loc->getMapForBlock(1)->getLocalNumElements());
for (size_type j = 0; j < data.size(); j++)
{
double target = .456 + myRank + i;
TEST_ASSERT(data[j] == target || data[j] == 2.0 * target);
}
Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorBase<double>>(param_f_vec->getVectorBlock(2))->getLocalData(Teuchos::ptrFromRef(data));
TEST_EQUALITY(static_cast<size_t>(data.size()), b_loc->getMapForBlock(2)->getLocalNumElements());
for (size_type j = 0; j < data.size(); j++)
{
double target = .789 + myRank + i;
TEST_ASSERT(data[j] == target || data[j] == 2.0 * target);
}
}
}

Teuchos::RCP<panzer::PureBasis> buildBasis(std::size_t worksetSize, const std::string &basisName)
{
Teuchos::RCP<shards::CellTopology> topo =
Expand Down
Loading

0 comments on commit 49a91f1

Please sign in to comment.