Panzer  Version of the Day
Panzer_Integrator_BasisTimesVector_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef __Panzer_Integerator_BasisTimesVector_impl_hpp__
44 #define __Panzer_Integerator_BasisTimesVector_impl_hpp__
45 
46 #include "Intrepid2_FunctionSpaceTools.hpp"
48 #include "Panzer_BasisIRLayout.hpp"
51 #include "Phalanx_KokkosDeviceTypes.hpp"
52 
53 namespace panzer {
54 
55 //**********************************************************************
56 PHX_EVALUATOR_CTOR(Integrator_BasisTimesVector,p) :
57  residual( p.get<std::string>("Residual Name"),
58  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
59  vectorField( p.get<std::string>("Value Name"),
60  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector),
61  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
62 {
64  p.validateParameters(*valid_params);
65 
67  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
68 
69  // Verify that this basis supports the gradient operation
70  TEUCHOS_TEST_FOR_EXCEPTION(!basis->isVectorBasis(),std::logic_error,
71  "Integrator_BasisTimesVector: Basis of type \"" << basis->name() << "\" is not "
72  "a vector basis.");
74  "Integrator_BasisTimesVector: Basis of type \"" << basis->name() << "\" does not "
75  "require orientation. This seems very strange, so I'm failing.");
76 
77  this->addEvaluatedField(residual);
78  this->addDependentField(vectorField);
79 
80  multiplier = p.get<double>("Multiplier");
81 
82 
83  if (p.isType<Teuchos::RCP<const std::vector<std::string> > >("Field Multipliers")) {
84  const std::vector<std::string>& field_multiplier_names =
85  *(p.get<Teuchos::RCP<const std::vector<std::string> > >("Field Multipliers"));
86 
87  int i=0;
88  field_multipliers.resize(field_multiplier_names.size());
89  kokkos_field_multipliers = Kokkos::View<Kokkos::View<ScalarT**>* >("BasisTimesVector::KokkosFieldMultipliers", field_multiplier_names.size());
90  Kokkos::fence();
91  for (std::vector<std::string>::const_iterator name =
92  field_multiplier_names.begin();
93  name != field_multiplier_names.end(); ++name) {
94  PHX::MDField<ScalarT,Cell,IP> tmp_field(*name, p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
95  Kokkos::fence();
96  field_multipliers[i++] = tmp_field;
97  this->addDependentField(field_multipliers[i-1]);
98  }
99  }
100 
101  std::string n = "Integrator_BasisTimesVector: " + residual.fieldTag().name();
102  this->setName(n);
103 }
104 
105 //**********************************************************************
106 PHX_POST_REGISTRATION_SETUP(Integrator_BasisTimesVector,sd,fm)
107 {
108  this->utils.setFieldData(residual,fm);
109  this->utils.setFieldData(vectorField,fm);
110  // this->utils.setFieldData(dof_orientation,fm);
111 
112  for (std::size_t i=0; i<field_multipliers.size(); ++i) {
113  this->utils.setFieldData(field_multipliers[i],fm);
114  kokkos_field_multipliers(i) = field_multipliers[i].get_static_view();
115  }
116 
117  basis_card = residual.dimension(1); // basis cardinality
118  num_qp = vectorField.dimension(1);
119  num_dim = vectorField.dimension(2); // dimension of a vector
120 
121  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
122 
123  // tmp = Intrepid2:FieldContainer<ScalarT>(vectorField.dimension(0), num_qp, num_dim);
124  MDFieldArrayFactory af("",fm.template getKokkosExtendedDataTypeDimensions<EvalT>(),true);
125 }
126 
127 
128 //**********************************************************************
129 template<typename EvalT, typename TRAITS>
130 template<int NUM_FIELD_MULT>
131 KOKKOS_INLINE_FUNCTION
132 void Integrator_BasisTimesVector<EvalT, TRAITS>::operator()(const FieldMultTag<NUM_FIELD_MULT> &, const size_t &cell) const {
133  const int nqp = vectorField.extent_int(1), ndim = vectorField.extent_int(2);
134  const int nfm = kokkos_field_multipliers.extent_int(0);
135  const int nbf = weighted_basis_vector.extent_int(1);
136 
137  for (int lbf = 0; lbf < nbf; lbf++)
138  residual(cell,lbf) = 0.0;
139 
140  ScalarT tmp, fmm=1;
141  if ( NUM_FIELD_MULT == 0 ){
142  for (int qp = 0; qp < nqp; ++qp) {
143  for (int d = 0; d < ndim; ++d) {
144  tmp = multiplier * vectorField(cell,qp,d);
145  for (int lbf = 0; lbf < nbf; lbf++)
146  residual(cell,lbf) += weighted_basis_vector(cell, lbf, qp, d)*tmp;
147  }
148  }
149  } else if ( NUM_FIELD_MULT == 1 ){
150  for (int qp = 0; qp < nqp; ++qp) {
151  for (int d = 0; d < ndim; ++d) {
152  tmp = multiplier * vectorField(cell,qp,d)*kokkos_field_multipliers(0)(cell,qp);
153  for (int lbf = 0; lbf < nbf; lbf++)
154  residual(cell,lbf) += weighted_basis_vector(cell, lbf, qp, d)*tmp;
155  }
156  }
157  } else {
158  for (int qp = 0; qp < nqp; ++qp) {
159  for (int i = 0; i < nfm; ++i)
160  fmm *= kokkos_field_multipliers(i)(cell,qp);
161  for (int d = 0; d < ndim; ++d) {
162  tmp = multiplier * vectorField(cell,qp,d)*fmm;
163  for (int lbf = 0; lbf < nbf; lbf++)
164  residual(cell,lbf) += weighted_basis_vector(cell, lbf, qp, d)*tmp;
165  }
166  }
167  }
168 }
169 
170 //**********************************************************************
171 PHX_EVALUATE_FIELDS(Integrator_BasisTimesVector,workset)
172 {
173  // residual.deep_copy(ScalarT(0.0));
174  weighted_basis_vector = this->wda(workset).bases[basis_index]->weighted_basis_vector;
175  if ( field_multipliers.size() == 0)
176  Kokkos::parallel_for(Kokkos::RangePolicy<FieldMultTag<0> >(0, workset.num_cells),*this);
177  else if ( field_multipliers.size() == 1)
178  Kokkos::parallel_for(Kokkos::RangePolicy<FieldMultTag<1> >(0, workset.num_cells),*this);
179  else
180  Kokkos::parallel_for(Kokkos::RangePolicy<FieldMultTag<-1> >(0, workset.num_cells),*this);
181 }
182 
183 //**********************************************************************
184 template<typename EvalT, typename TRAITS>
187 {
189  p->set<std::string>("Residual Name", "?");
190  p->set<std::string>("Value Name", "?");
191  p->set<std::string>("Test Field Name", "?");
193  p->set("Basis", basis);
195  p->set("IR", ir);
196  p->set<double>("Multiplier", 1.0);
198  p->set("Field Multipliers", fms);
199  return p;
200 }
201 
202 //**********************************************************************
203 
204 }
205 
206 #endif
207 
PHX::MDField< ScalarT > residual
Evaluates a Dirichlet BC residual corresponding to a field value.
std::size_t basis_index
std::vector< std::string >::size_type getBasisIndex(std::string basis_name, panzer::Workset &workset, WorksetDetailsAccessor &wda)
Returns the index in the workset bases for a particular BasisIRLayout name.
std::string name() const
A unique key that is the combination of the basis type and basis order.
std::string basis_name
PHX::MDField< ScalarT, Cell > tmp
std::size_t num_qp
Kokkos::View< Kokkos::View< ScalarT **> *> kokkos_field_multipliers
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
PHX::MDField< ScalarT, Cell, IP, Dim > vectorField
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void operator()(const FieldMultTag< NUM_FIELD_MULT > &, const std::size_t &cell) const
T * get() const
bool isVectorBasis() const
bool requiresOrientations() const
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
PHX::MDField< double, Cell, BASIS, IP, Dim > weighted_basis_vector
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
PHX_EVALUATOR_CTOR(BasisValues_Evaluator, p)
double multiplier
PHX_EVALUATE_FIELDS(BasisValues_Evaluator, workset)
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
std::vector< PHX::MDField< const ScalarT, Cell, IP > > field_multipliers
PHX_POST_REGISTRATION_SETUP(BasisValues_Evaluator, sd, fm)