Albany: a Trilinos-based PDE code: PHAL_NSMaterialProperty

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <fstream>
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 #include "Sacado_ParameterRegistration.hpp"
00011 #include "Albany_Utils.hpp"
00012 
00013 namespace PHAL {
00014 
00015 template<typename EvalT, typename Traits>
00016 NSMaterialProperty<EvalT, Traits>::
00017 NSMaterialProperty(Teuchos::ParameterList& p) :
00018   name_mp(p.get<std::string>("Material Property Name")),
00019   layout(p.get<Teuchos::RCP<PHX::DataLayout> >("Data Layout")),
00020   matprop(name_mp,layout),
00021   rank(layout->rank()),
00022   dims(),
00023   matPropType(SCALAR_CONSTANT)
00024 {
00025   layout->dimensions(dims);
00026 
00027   double default_value = p.get("Default Value", 1.0);
00028 
00029   Teuchos::RCP<ParamLib> paramLib = 
00030     p.get< Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00031 
00032   Teuchos::ParameterList* mp_list = 
00033     p.get<Teuchos::ParameterList*>("Parameter List");
00034   std::string type = mp_list->get("Type", "Constant");
00035   if (type == "Constant") {
00036     if (rank == 2) {
00037       matPropType = SCALAR_CONSTANT;
00038       scalar_constant_value = mp_list->get("Value", default_value);
00039 
00040       // Add property as a Sacado-ized parameter
00041       new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00042   name_mp, this, paramLib);
00043     }
00044     else if (rank == 3) {
00045       matPropType = VECTOR_CONSTANT;
00046       PHX::DataLayout::size_type numDims = dims[2];
00047       Teuchos::Array<double> tmp = 
00048   mp_list->get< Teuchos::Array<double> >("Value");
00049       vector_constant_value.resize(numDims);
00050       TEUCHOS_TEST_FOR_EXCEPTION(vector_constant_value.size() != numDims,
00051        std::logic_error,
00052        "Vector constant value for material property " <<
00053        name_mp << " has size " << 
00054        vector_constant_value.size() << " but expected size "
00055        << numDims);
00056 
00057       for (PHX::DataLayout::size_type i=0; i<numDims; i++)
00058   vector_constant_value[i] = tmp[i];
00059 
00060       // Add property as a Sacado-ized parameter
00061       for (PHX::DataLayout::size_type i=0; i<numDims; i++)
00062   new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00063     Albany::strint(name_mp,i), this, paramLib);
00064     }
00065     else if (rank == 4) {
00066       matPropType = TENSOR_CONSTANT;
00067       PHX::DataLayout::size_type numRows = dims[2];
00068       PHX::DataLayout::size_type numCols = dims[3];
00069       Teuchos::TwoDArray<double> tmp = 
00070   mp_list->get< Teuchos::TwoDArray<double> >("Value");
00071       TEUCHOS_TEST_FOR_EXCEPTION(tensor_constant_value.getNumRows() != numRows ||
00072        tensor_constant_value.getNumCols() != numCols,
00073        std::logic_error,
00074        "Tensor constant value for material property " <<
00075        name_mp << " has dimensions " << 
00076        tensor_constant_value.getNumRows() << "x" <<
00077        tensor_constant_value.getNumCols() << 
00078        " but expected dimensions " << 
00079        numRows << "x" << numCols);
00080       tensor_constant_value = Teuchos::TwoDArray<ScalarT>(numRows, numCols);
00081       for (PHX::DataLayout::size_type i=0; i<numRows; i++)
00082   for (PHX::DataLayout::size_type j=0; j<numCols; j++)
00083     tensor_constant_value(i,j) = tmp(i,j);
00084 
00085       // Add property as a Sacado-ized parameter
00086       for (PHX::DataLayout::size_type i=0; i<numRows; i++)
00087   for (PHX::DataLayout::size_type j=0; j<numCols; j++)
00088     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00089       Albany::strint(Albany::strint(name_mp,i),j), this, paramLib);
00090     }
00091     else
00092       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00093        "Invalid material property rank " << rank << 
00094        ".  Acceptable values are 2 (scalar), " << 
00095        "3 (vector), or 4 (tensor)");
00096   }
00097   else if (type == "Truncated KL Expansion" || 
00098      type == "Log Normal RF" ||
00099      type == "Exponential Truncated KL Expansion") {
00100     if (type == "Truncated KL Expansion")
00101       matPropType = KL_RAND_FIELD;
00102     else if (type == "Log Normal RF" || 
00103        type == "Exponential Truncated KL Expansion")
00104       matPropType = EXP_KL_RAND_FIELD;
00105  
00106     Teuchos::RCP<PHX::DataLayout> coord_dl = 
00107       p.get< Teuchos::RCP<PHX::DataLayout> >("Coordinate Vector Data Layout");
00108     coordVec = PHX::MDField<MeshScalarT>(
00109       p.get<std::string>("Coordinate Vector Name"),
00110       coord_dl);
00111     this->addDependentField(coordVec);
00112     std::vector<PHX::DataLayout::size_type> coord_dims;
00113     coord_dl->dimensions(coord_dims);
00114     point.resize(coord_dims[2]);
00115 
00116     exp_rf_kl = 
00117       Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<MeshScalarT>(*mp_list));
00118     int num_KL = exp_rf_kl->stochasticDimension();
00119 
00120     // Add KL random variables as Sacado-ized parameters
00121     rv.resize(num_KL);
00122     for (int i=0; i<num_KL; i++) {
00123       std::string ss = Albany::strint(name_mp + " KL Random Variable",i);
00124       new Sacado::ParameterRegistration<EvalT, SPL_Traits>(ss, this, paramLib);
00125       rv[i] = mp_list->get(ss, 0.0);
00126     }
00127   }
00128   else if (type == "SQRT Temperature Dependent") {
00129     matPropType = SQRT_TEMP;
00130     scalar_constant_value = mp_list->get("Reference Value", default_value);
00131     ref_temp = mp_list->get("Reference Temperature", default_value);
00132     T = PHX::MDField<ScalarT>(
00133       p.get<std::string>("Temperature Variable Name"),
00134       layout);
00135     this->addDependentField(T); 
00136 
00137     // Add property as a Sacado-ized parameter
00138     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00139       name_mp+" Reference Value", this, paramLib);
00140   }
00141   else if (type == "invSQRT Temperature Dependent") {
00142     matPropType = INV_SQRT_TEMP;
00143     scalar_constant_value = mp_list->get("Reference Value", default_value);
00144     ref_temp = mp_list->get("Reference Temperature", default_value);
00145     T = PHX::MDField<ScalarT>(
00146       p.get<std::string>("Temperature Variable Name"),
00147       layout);
00148     this->addDependentField(T); 
00149 
00150     // Add property as a Sacado-ized parameter
00151     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00152       name_mp+" Reference Value", this, paramLib);
00153   }
00154   else if (type == "Transport Mean Free Path") {
00155     matPropType = NEUTRON_DIFFUSION;
00156     sigma_a = PHX::MDField<ScalarT>(
00157       p.get<std::string>("Absorption Cross Section Name"),
00158       layout);
00159     sigma_s = PHX::MDField<ScalarT>(
00160       p.get<std::string>("Scattering Cross Section Name"),
00161       layout);
00162     mu = PHX::MDField<ScalarT>(
00163       p.get<std::string>("Average Scattering Angle Name"),
00164       layout);
00165     this->addDependentField(sigma_a); 
00166     this->addDependentField(sigma_s); 
00167     this->addDependentField(mu);
00168   }
00169   else {
00170     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00171            "Invalid material property type " << type);
00172   } 
00173 
00174   this->addEvaluatedField(matprop);
00175   this->setName(name_mp+PHX::TypeString<EvalT>::value);
00176 }
00177 
00178 // **********************************************************************
00179 template<typename EvalT, typename Traits>
00180 void NSMaterialProperty<EvalT, Traits>::
00181 postRegistrationSetup(typename Traits::SetupData d,
00182                       PHX::FieldManager<Traits>& fm)
00183 {
00184   this->utils.setFieldData(matprop,fm);
00185   if (matPropType == KL_RAND_FIELD || matPropType == EXP_KL_RAND_FIELD) 
00186     this->utils.setFieldData(coordVec,fm);
00187   if (matPropType == SQRT_TEMP || matPropType == INV_SQRT_TEMP) 
00188     this->utils.setFieldData(T,fm);
00189   if (matPropType == NEUTRON_DIFFUSION) {
00190     this->utils.setFieldData(sigma_a,fm);
00191     this->utils.setFieldData(sigma_s,fm);
00192     this->utils.setFieldData(mu,fm);
00193   }
00194 }
00195 
00196 // **********************************************************************
00197 template<typename EvalT, typename Traits>
00198 void NSMaterialProperty<EvalT, Traits>::
00199 evaluateFields(typename Traits::EvalData workset)
00200 {
00201   if (matPropType == SCALAR_CONSTANT) {
00202     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00203       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00204   matprop(cell,qp) = scalar_constant_value;
00205       }
00206     }
00207   }
00208   else if (matPropType == VECTOR_CONSTANT) {
00209     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00210       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00211   for (std::size_t dim=0; dim < dims[2]; ++dim) {
00212     matprop(cell,qp,dim) = vector_constant_value[dim];
00213   }
00214       }
00215     }
00216   }
00217   else if (matPropType == TENSOR_CONSTANT) {
00218     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00219       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00220   for (std::size_t dim1=0; dim1 < dims[2]; ++dim1) {
00221     for (std::size_t dim2=0; dim2 < dims[3]; ++dim2) {
00222       matprop(cell,qp,dim1,dim2) = tensor_constant_value(dim1,dim2);
00223     }
00224   }
00225       }
00226     }
00227   }
00228   else if (matPropType == SQRT_TEMP) {
00229     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00230       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00231   matprop(cell,qp) = scalar_constant_value / sqrt(ref_temp) * sqrt(T(cell,qp));
00232       } 
00233     }
00234   }
00235   else if (matPropType == INV_SQRT_TEMP) {
00236     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00237       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00238   matprop(cell,qp) = scalar_constant_value * sqrt(ref_temp) / sqrt(T(cell,qp));
00239       } 
00240     }
00241   }
00242   else if (matPropType == NEUTRON_DIFFUSION) {
00243     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00244       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00245   matprop(cell,qp) = 
00246     1.0 / (3.0*(1.0 - mu(cell,qp))*(sigma_a(cell,qp) + sigma_s(cell,qp)));
00247       } 
00248     }
00249   }
00250   else {
00251     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00252       for (std::size_t qp=0; qp < dims[1]; ++qp) {
00253   for (std::size_t i=0; i<point.size(); i++)
00254     point[i] = 
00255       Sacado::ScalarValue<MeshScalarT>::eval(coordVec(cell,qp,i));
00256   matprop(cell,qp) = exp_rf_kl->evaluate(point, rv);       
00257         if (matPropType == EXP_KL_RAND_FIELD)
00258           matprop(cell,qp) = std::exp(matprop(cell,qp));       
00259       }
00260     }
00261   }
00262 }
00263 
00264 // **********************************************************************
00265 template<typename EvalT,typename Traits>
00266 typename NSMaterialProperty<EvalT,Traits>::ScalarT& 
00267 NSMaterialProperty<EvalT,Traits>::getValue(const std::string &n)
00268 {
00269   if (matPropType == SCALAR_CONSTANT || 
00270       matPropType == SQRT_TEMP || 
00271       matPropType == INV_SQRT_TEMP) {
00272     return scalar_constant_value;
00273   }
00274   else if (matPropType == VECTOR_CONSTANT) {
00275     for (std::size_t dim=0; dim<vector_constant_value.size(); ++dim)
00276       if (n == Albany::strint(name_mp,dim))
00277   return vector_constant_value[dim];
00278   }
00279   else if (matPropType == TENSOR_CONSTANT) {
00280     for (std::size_t dim1=0; dim1<tensor_constant_value.getNumRows(); ++dim1)
00281       for (std::size_t dim2=0; dim2<tensor_constant_value.getNumCols(); ++dim2)
00282   if (n == Albany::strint(Albany::strint(name_mp,dim1),dim2))
00283     return tensor_constant_value(dim1,dim2);
00284   }
00285   else if (matPropType == KL_RAND_FIELD || matPropType == EXP_KL_RAND_FIELD) {
00286     for (int i=0; i<rv.size(); i++)
00287       if (n == Albany::strint(name_mp + " KL Random Variable",i))
00288   return rv[i];
00289   }
00290   TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00291          std::endl <<
00292          "Error! Logic error in getting paramter " << n
00293          << " in NSMaterialProperty::getValue()" << std::endl);
00294   return scalar_constant_value;
00295 }
00296 
00297 // **********************************************************************
00298 // **********************************************************************
00299 }
00300
PHAL_NSMaterialProperty_Def.hpp
