Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Panzer tangent scatter blocked tpetra #13601

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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
Loading