Skip to content

Commit

Permalink
Panzer: Lazy evaluate integration values 2 from workset (#9709)
Browse files Browse the repository at this point in the history 
Adding a lazy evaluation construction path to IntegrationValues2
  • Loading branch information
@seamill
seamill authored Sep 28, 2021
1 parent ce1a604 commit 7527310
Show file tree
Hide file tree
Showing 14 changed files with 2,066 additions and 1,494 deletions.
37 changes: 20 additions & 17 deletions packages/panzer/adapters-stk/example/PoissonInterfaceExample/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,7 @@
#include "Panzer_BlockedEpetraLinearObjFactory.hpp"
#include "Panzer_DOFManagerFactory.hpp"
#include "Panzer_DOFManager.hpp"
#include "Panzer_OrientationsInterface.hpp"
#include "Panzer_FieldManagerBuilder.hpp"
#include "Panzer_PureBasis.hpp"
#include "Panzer_GlobalData.hpp"
Expand Down Expand Up @@ -568,23 +569,6 @@ int main (int argc, char* argv[])
mesh_factory->completeMeshConstruction(*mesh,MPI_COMM_WORLD);
}

// build worksets
////////////////////////////////////////////////////////

Teuchos::RCP<panzer_stk::WorksetFactory> wkstFactory
= Teuchos::rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory
Teuchos::RCP<panzer::WorksetContainer> wkstContainer // attach it to a workset container (uses lazy evaluation)
= Teuchos::rcp(new panzer::WorksetContainer);
wkstContainer->setFactory(wkstFactory);
for(size_t i=0;i<physicsBlocks.size();i++)
wkstContainer->setNeeds(physicsBlocks[i]->elementBlockID(),physicsBlocks[i]->getWorksetNeeds());
wkstContainer->setWorksetSize(workset_size);

std::vector<std::string> elementBlockNames;
mesh->getElementBlockNames(elementBlockNames);
std::map<std::string,Teuchos::RCP<std::vector<panzer::Workset> > > volume_worksets;
panzer::getVolumeWorksetsFromContainer(*wkstContainer,elementBlockNames,volume_worksets);

// build DOF Manager and linear object factory
/////////////////////////////////////////////////////////////

Expand All @@ -603,6 +587,25 @@ int main (int argc, char* argv[])
checkInterfaceConnections(conn_manager, dofManager->getComm());
}

// build worksets
////////////////////////////////////////////////////////

Teuchos::RCP<panzer_stk::WorksetFactory> wkstFactory
= Teuchos::rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory
wkstFactory->setOrientationsInterface(Teuchos::rcp(new panzer::OrientationsInterface(dofManager)));
Teuchos::RCP<panzer::WorksetContainer> wkstContainer // attach it to a workset container (uses lazy evaluation)
= Teuchos::rcp(new panzer::WorksetContainer);
wkstContainer->setFactory(wkstFactory);
for(size_t i=0;i<physicsBlocks.size();i++)
wkstContainer->setNeeds(physicsBlocks[i]->elementBlockID(),physicsBlocks[i]->getWorksetNeeds());
wkstContainer->setWorksetSize(workset_size);

std::vector<std::string> elementBlockNames;
mesh->getElementBlockNames(elementBlockNames);
std::map<std::string,Teuchos::RCP<std::vector<panzer::Workset> > > volume_worksets;
panzer::getVolumeWorksetsFromContainer(*wkstContainer,elementBlockNames,volume_worksets);


// construct some linear algebra object, build object to pass to evaluators
Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > linObjFactory
= Teuchos::rcp(new panzer::BlockedEpetraLinearObjFactory<panzer::Traits,int>(tComm.getConst(),dofManager));
Expand Down
8 changes: 3 additions & 5 deletions packages/panzer/adapters-stk/test/panzer_workset_builder/check_rotation_matrices.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,13 +120,11 @@ namespace panzer {
//////////////////////////////////////////////////////////////

panzer::IntegrationDescriptor sid(2*2, panzer::IntegrationDescriptor::SURFACE);
std::map<std::string, panzer::WorksetNeeds> wkstRequirements;
wkstRequirements[element_block].addIntegrator(sid);

RCP<panzer_stk::WorksetFactory> wkstFactory
= rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory
RCP<panzer::WorksetContainer> wkstContainer // attach it to a workset container (uses lazy evaluation)
= rcp(new panzer::WorksetContainer(wkstFactory,wkstRequirements));
= rcp(new panzer::WorksetContainer(wkstFactory));

wkstContainer->setGlobalIndexer(dof_manager);

Expand All @@ -138,8 +136,8 @@ namespace panzer {

auto & workset = (*worksets)[0];

auto rot_matrices = workset.getIntegrationValues(sid).surface_rotation_matrices;
auto normals = workset.getIntegrationValues(sid).surface_normals;
auto rot_matrices = workset.getIntegrationValues(sid).getSurfaceRotationMatrices();
auto normals = workset.getIntegrationValues(sid).getSurfaceNormals();

const int num_owned_cells = workset.numOwnedCells();
const int num_ghost_cells = workset.numGhostCells();
Expand Down
32 changes: 16 additions & 16 deletions packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,18 +271,18 @@ namespace panzer {
bool is_uniform_;

// Only valid for uniform reference space
PHX::MDField<Scalar,IP,Dim> cubature_points_uniform_ref_;
PHX::MDField<const Scalar,IP,Dim> cubature_points_uniform_ref_;

// For non-uniform reference space
PHX::MDField<Scalar,Cell,IP,Dim> cubature_points_ref_;
PHX::MDField<const Scalar,Cell,IP,Dim> cubature_points_ref_;

// Geometry objects that represent cubature/point values
PHX::MDField<Scalar,Cell,IP,Dim,Dim> cubature_jacobian_;
PHX::MDField<Scalar,Cell,IP> cubature_jacobian_determinant_;
PHX::MDField<Scalar,Cell,IP,Dim,Dim> cubature_jacobian_inverse_;
PHX::MDField<Scalar,Cell,IP> cubature_weights_;
PHX::MDField<const Scalar,Cell,IP,Dim,Dim> cubature_jacobian_;
PHX::MDField<const Scalar,Cell,IP> cubature_jacobian_determinant_;
PHX::MDField<const Scalar,Cell,IP,Dim,Dim> cubature_jacobian_inverse_;
PHX::MDField<const Scalar,Cell,IP> cubature_weights_;

PHX::MDField<Scalar,Cell,NODE,Dim> cell_vertex_coordinates_;
PHX::MDField<const Scalar,Cell,NODE,Dim> cell_vertex_coordinates_;

// Number of cells to apply orientations to (required in situations where virtual cells exist)
int num_orientations_cells_;
Expand Down Expand Up @@ -329,10 +329,10 @@ namespace panzer {
*/
void
setup(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
PHX::MDField<Scalar, Cell, IP, Dim> reference_points,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
PHX::MDField<const Scalar, Cell, IP, Dim> reference_points,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<const Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
const int num_evaluated_cells = -1);

/**
Expand All @@ -350,10 +350,10 @@ namespace panzer {
*/
void
setupUniform(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
PHX::MDField<Scalar, IP, Dim> reference_points,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
PHX::MDField<const Scalar, IP, Dim> reference_points,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<const Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
const int num_evaluated_cells = -1);

/// Set the orientations object for applying orientations using the lazy evaluation path - required for certain bases
Expand All @@ -363,7 +363,7 @@ namespace panzer {

/// Set the cubature weights (weighted measure) for the basis values object - required to get weighted basis objects
void
setWeightedMeasure(PHX::MDField<Scalar, Cell, IP> weighted_measure);
setWeightedMeasure(PHX::MDField<const Scalar, Cell, IP> weighted_measure);

/// Set the cell vertex coordinates (required for getBasisCoordinates())
void
Expand Down
78 changes: 52 additions & 26 deletions packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@
#include "Intrepid2_Orientation.hpp"
#include "Intrepid2_OrientationTools.hpp"

// FIXME: There are some calls in Intrepid2 that require non-const arrays when they should be const - search for PHX::getNonConstDynRankViewFromConstMDField
#include "Phalanx_GetNonConstDynRankViewFromConstMDField.hpp"

namespace panzer {
namespace {

Expand Down Expand Up @@ -599,10 +602,10 @@ template <typename Scalar>
void
BasisValues2<Scalar>::
setup(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
PHX::MDField<Scalar, Cell, IP, Dim> reference_points,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
PHX::MDField<const Scalar, Cell, IP, Dim> reference_points,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<const Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
const int num_evaluated_cells)
{
basis_layout = basis;
Expand All @@ -626,10 +629,10 @@ template <typename Scalar>
void
BasisValues2<Scalar>::
setupUniform(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
PHX::MDField<Scalar, IP, Dim> reference_points,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
PHX::MDField<const Scalar, IP, Dim> reference_points,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian,
PHX::MDField<const Scalar, Cell, IP> point_jacobian_determinant,
PHX::MDField<const Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
const int num_evaluated_cells)
{
basis_layout = basis;
Expand Down Expand Up @@ -735,7 +738,7 @@ resetArrays()
template <typename Scalar>
void
BasisValues2<Scalar>::
setWeightedMeasure(PHX::MDField<Scalar, Cell, IP> weighted_measure)
setWeightedMeasure(PHX::MDField<const Scalar, Cell, IP> weighted_measure)
{
TEUCHOS_TEST_FOR_EXCEPT_MSG(build_weighted,
"BasisValues2::setWeightedMeasure : Weighted measure already set. Can only set weighted measure once after setup or setupUniform have beens called.");
Expand Down Expand Up @@ -807,7 +810,9 @@ getBasisValuesRef(const bool cache,
const int num_card = basis_layout->cardinality();

auto tmp_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_quad);
intrepid_basis->getValues(tmp_basis_ref_scalar.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_VALUE);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(tmp_basis_ref_scalar.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_VALUE);
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -840,7 +845,9 @@ getVectorBasisValuesRef(const bool cache,
const int num_dim = basis_layout->dimension();

auto tmp_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
intrepid_basis->getValues(tmp_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(tmp_basis_ref_vector.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_VALUE);
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -873,7 +880,9 @@ getGradBasisValuesRef(const bool cache,
const int num_dim = basis_layout->dimension();

auto tmp_grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
intrepid_basis->getValues(tmp_grad_basis_ref.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_GRAD);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(tmp_grad_basis_ref.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_GRAD);
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -906,7 +915,9 @@ getCurl2DVectorBasisRef(const bool cache,
const int num_card = basis_layout->cardinality();

auto tmp_curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("dyn_curl_basis_ref_scalar",num_card,num_quad);
intrepid_basis->getValues(tmp_curl_basis_ref_scalar.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_CURL);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(tmp_curl_basis_ref_scalar.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_CURL);
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -940,7 +951,9 @@ getCurlVectorBasisRef(const bool cache,
const int num_dim = basis_layout->dimension();

auto tmp_curl_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim);
intrepid_basis->getValues(tmp_curl_basis_ref_vector.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_CURL);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(tmp_curl_basis_ref_vector.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_CURL);
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -972,7 +985,9 @@ getDivVectorBasisRef(const bool cache,
const int num_card = basis_layout->cardinality();

auto tmp_div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("dyn_div_basis_ref_scalar",num_card,num_quad);
intrepid_basis->getValues(tmp_div_basis_ref.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_DIV);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(tmp_div_basis_ref.get_view(), cubature_points_uniform_ref, Intrepid2::OPERATOR_DIV);
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -1004,13 +1019,12 @@ getBasisCoordinates(const bool cache,
auto tmp_basis_coordinates = af.buildStaticArray<Scalar, Cell, BASIS, IP>("basis_coordinates", num_cells_, num_card, num_dim);
auto s_aux = Kokkos::subview(tmp_basis_coordinates.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());

// Some kind of type bug in Intrepid's mapToPhysicalFrame call? Can't let bcr be const
auto bcr = getBasisCoordinatesRef(false);
Kokkos::DynRankView<double> bcr_copy("nonconst_bcr",bcr.extent(0),bcr.extent(1));
Kokkos::deep_copy(bcr_copy,bcr.get_view());
// Don't forget that since we are not caching this, we have to make sure the managed view remains alive while we use the non-const wrapper
auto const_bcr = getBasisCoordinatesRef(false);
auto bcr = PHX::getNonConstDynRankViewFromConstMDField(const_bcr);

Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
cell_tools.mapToPhysicalFrame(s_aux, bcr_copy, s_vertex_coordinates, intrepid_basis->getBaseCellTopology());
cell_tools.mapToPhysicalFrame(s_aux, bcr, s_vertex_coordinates, intrepid_basis->getBaseCellTopology());
PHX::Device().fence();

// Store for later if cache is enabled
Expand Down Expand Up @@ -1081,8 +1095,10 @@ getBasisValues(const bool weighted,

if(hasUniformReferenceSpace()){

auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

// Apply a single reference representation to all cells
intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE);
intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_VALUE);

const std::pair<int,int> cell_range(0,num_evaluate_cells_);
auto s_aux = Kokkos::subview(tmp_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
Expand Down Expand Up @@ -1239,8 +1255,10 @@ getVectorBasisValues(const bool weighted,

if(hasUniformReferenceSpace()){

auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

// Apply a single reference representation to all cells
intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE);
intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_VALUE);

const std::pair<int,int> cell_range(0,num_evaluate_cells_);
auto s_aux = Kokkos::subview(tmp_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
Expand Down Expand Up @@ -1414,8 +1432,10 @@ getGradBasisValues(const bool weighted,

if(hasUniformReferenceSpace()){

auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

// Apply a single reference representation to all cells
intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_GRAD);
intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_GRAD);

const std::pair<int,int> cell_range(0,num_evaluate_cells_);
auto s_aux = Kokkos::subview(tmp_grad_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
Expand Down Expand Up @@ -1560,7 +1580,9 @@ getCurl2DVectorBasis(const bool weighted,

if(hasUniformReferenceSpace()){

intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_CURL);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_CURL);

const std::pair<int,int> cell_range(0,num_evaluate_cells_);
auto s_aux = Kokkos::subview(tmp_curl_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
Expand Down Expand Up @@ -1703,7 +1725,9 @@ getCurlVectorBasis(const bool weighted,

if(hasUniformReferenceSpace()){

intrepid_basis->getValues(cell_curl_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_CURL);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(cell_curl_basis_ref_vector.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_CURL);

const std::pair<int,int> cell_range(0,num_evaluate_cells_);
auto s_aux = Kokkos::subview(tmp_curl_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
Expand Down Expand Up @@ -1848,7 +1872,9 @@ getDivVectorBasis(const bool weighted,

if(hasUniformReferenceSpace()){

intrepid_basis->getValues(cell_div_basis_ref.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_DIV);
auto cubature_points_uniform_ref = PHX::getNonConstDynRankViewFromConstMDField(cubature_points_uniform_ref_);

intrepid_basis->getValues(cell_div_basis_ref.get_view(),cubature_points_uniform_ref,Intrepid2::OPERATOR_DIV);

const std::pair<int,int> cell_range(0,num_evaluate_cells_);
auto s_aux = Kokkos::subview(tmp_div_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
Expand Down
2 changes: 1 addition & 1 deletion packages/panzer/disc-fe/src/Panzer_IntegrationDescriptor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ IntegrationDescriptor::setup(const int cubature_order, const int integration_typ
_cubature_order = cubature_order;
_side = side;

if(_integration_type == SIDE or _integration_type == CV_BOUNDARY){
if((_integration_type == SIDE) or (_integration_type == CV_BOUNDARY)){
TEUCHOS_ASSERT(side >= 0);
} else {
TEUCHOS_ASSERT(side == -1);
Expand Down
Loading

0 comments on commit 7527310

Please sign in to comment.