Albany: a Trilinos-based PDE code: PHAL_ThermalConductivity

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 //Radom field types
00014 enum SG_RF {CONSTANT, UNIFORM, LOGNORMAL};
00015 const int num_sg_rf = 3;
00016 const SG_RF sg_rf_values[] = {CONSTANT, UNIFORM, LOGNORMAL};
00017 const char *sg_rf_names[] = {"Constant", "Uniform", "Log-Normal"};
00018 
00019 SG_RF randField = CONSTANT;
00020 
00021 namespace PHAL {
00022 
00023 template<typename EvalT, typename Traits>
00024 ThermalConductivity<EvalT, Traits>::
00025 ThermalConductivity(Teuchos::ParameterList& p) :
00026   thermalCond(p.get<std::string>("QP Variable Name"),
00027         p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"))
00028 {
00029 
00030   Teuchos::ParameterList* cond_list = 
00031     p.get<Teuchos::ParameterList*>("Parameter List");
00032 
00033   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00034                this->getValidThermalCondParameters();
00035 
00036   // Check the parameters contained in the input file. Do not check the defaults
00037   // set programmatically
00038   cond_list->validateParameters(*reflist, 0, 
00039     Teuchos::VALIDATE_USED_ENABLED, Teuchos::VALIDATE_DEFAULTS_DISABLED);
00040 
00041   Teuchos::RCP<PHX::DataLayout> vector_dl =
00042     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00043   std::vector<PHX::DataLayout::size_type> dims;
00044   vector_dl->dimensions(dims);
00045   numQPs  = dims[1];
00046   numDims = dims[2];
00047 
00048   std::string ebName = 
00049     p.get<std::string>("Element Block Name", "Missing");
00050 
00051   type = cond_list->get("Thermal Conductivity Type", "Constant");
00052 
00053   if (type == "Constant") {
00054 
00055     ScalarT value = cond_list->get("Value", 1.0);
00056     init_constant(value, p);
00057 
00058   } 
00059 
00060   else if (type == "Truncated KL Expansion" || type == "Log Normal RF") {
00061 
00062     init_KL_RF(type, *cond_list, p);
00063 
00064   } 
00065 
00066   else if (type == "Block Dependent") 
00067   {
00068     // We have a multiple material problem and need to map element blocks to material data
00069 
00070     if(p.isType<Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB")){
00071        materialDB = p.get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
00072     }
00073     else {
00074        TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00075          std::endl <<
00076          "Error! Must specify a material database if using block dependent " << 
00077          "thermal conductivity" << std::endl);
00078     }
00079 
00080     // Get the sublist for thermal conductivity for the element block in the mat DB (the material in the
00081     // elem block ebName.
00082 
00083     Teuchos::ParameterList& subList = materialDB->getElementBlockSublist(ebName, "Thermal Conductivity");
00084 
00085     std::string typ = subList.get("Thermal Conductivity Type", "Constant");
00086 
00087     if (typ == "Constant") {
00088 
00089        ScalarT value = subList.get("Value", 1.0);
00090        init_constant(value, p);
00091 
00092     }
00093     else if (typ == "Truncated KL Expansion" || typ == "Log Normal RF") {
00094 
00095        init_KL_RF(typ, subList, p);
00096 
00097     } 
00098   } // Block dependent
00099 
00100   else {
00101     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00102            "Invalid thermal conductivity type " << type);
00103   } 
00104 
00105   this->addEvaluatedField(thermalCond);
00106   this->setName("Thermal Conductivity"+PHX::TypeString<EvalT>::value);
00107 }
00108 
00109 template<typename EvalT, typename Traits>
00110 void
00111 ThermalConductivity<EvalT, Traits>::
00112 init_constant(ScalarT value, Teuchos::ParameterList& p){
00113 
00114     is_constant = true;
00115     randField = CONSTANT;
00116 
00117     constant_value = value;
00118 
00119     // Add thermal conductivity as a Sacado-ized parameter
00120     Teuchos::RCP<ParamLib> paramLib = 
00121       p.get< Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00122 
00123     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00124       "Thermal Conductivity", this, paramLib);
00125 
00126 } // init_constant
00127 
00128 template<typename EvalT, typename Traits>
00129 void
00130 ThermalConductivity<EvalT, Traits>::
00131 init_KL_RF(std::string &type, Teuchos::ParameterList& sublist, Teuchos::ParameterList& p){
00132 
00133     is_constant = false;
00134 
00135     if (type == "Truncated KL Expansion")
00136       randField = UNIFORM;
00137     else if (type == "Log Normal RF")
00138       randField = LOGNORMAL;
00139  
00140     Teuchos::RCP<PHX::DataLayout> scalar_dl =
00141       p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00142     Teuchos::RCP<PHX::DataLayout> vector_dl =
00143       p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00144     PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>
00145       fx(p.get<std::string>("QP Coordinate Vector Name"), vector_dl);
00146     coordVec = fx;
00147     this->addDependentField(coordVec);
00148 
00149     exp_rf_kl = 
00150       Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<MeshScalarT>(sublist));
00151     int num_KL = exp_rf_kl->stochasticDimension();
00152 
00153     // Add KL random variables as Sacado-ized parameters
00154     rv.resize(num_KL);
00155     Teuchos::RCP<ParamLib> paramLib = 
00156       p.get< Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00157     for (int i=0; i<num_KL; i++) {
00158       std::string ss = Albany::strint("Thermal Conductivity KL Random Variable",i);
00159       new Sacado::ParameterRegistration<EvalT, SPL_Traits>(ss, this, paramLib);
00160       rv[i] = sublist.get(ss, 0.0);
00161     }
00162 
00163 } // (type == "Truncated KL Expansion" || type == "Log Normal RF") 
00164 
00165 // **********************************************************************
00166 template<typename EvalT, typename Traits>
00167 void ThermalConductivity<EvalT, Traits>::
00168 postRegistrationSetup(typename Traits::SetupData d,
00169                       PHX::FieldManager<Traits>& fm)
00170 {
00171   this->utils.setFieldData(thermalCond,fm);
00172   if (!is_constant)
00173       this->utils.setFieldData(coordVec,fm);
00174 }
00175 
00176 // **********************************************************************
00177 template<typename EvalT, typename Traits>
00178 void ThermalConductivity<EvalT, Traits>::
00179 evaluateFields(typename Traits::EvalData workset)
00180 {
00181   if (is_constant) {
00182     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00183       for (std::size_t qp=0; qp < numQPs; ++qp) {
00184          thermalCond(cell,qp) = constant_value;
00185       }
00186     }
00187   }
00188 
00189   else {
00190     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00191       for (std::size_t qp=0; qp < numQPs; ++qp) {
00192           Teuchos::Array<MeshScalarT> point(numDims);
00193           for (std::size_t i=0; i<numDims; i++)
00194               point[i] = Sacado::ScalarValue<MeshScalarT>::eval(coordVec(cell,qp,i));
00195           if (randField == UNIFORM)
00196               thermalCond(cell,qp) = exp_rf_kl->evaluate(point, rv);       
00197           else if (randField == LOGNORMAL)
00198               thermalCond(cell,qp) = std::exp(exp_rf_kl->evaluate(point, rv));       
00199       }
00200     }
00201   }
00202 }
00203 
00204 // **********************************************************************
00205 template<typename EvalT,typename Traits>
00206 typename ThermalConductivity<EvalT,Traits>::ScalarT& 
00207 ThermalConductivity<EvalT,Traits>::getValue(const std::string &n)
00208 {
00209   if (is_constant) {
00210     return constant_value;
00211   }
00212 
00213   for (int i=0; i<rv.size(); i++) {
00214     if (n == Albany::strint("Thermal Conductivity KL Random Variable",i))
00215       return rv[i];
00216   }
00217   TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00218          std::endl <<
00219          "Error! Logic error in getting paramter " << n
00220          << " in ThermalConductivity::getValue()" << std::endl);
00221   return constant_value;
00222 }
00223 
00224 // **********************************************************************
00225 template<typename EvalT,typename Traits>
00226 Teuchos::RCP<const Teuchos::ParameterList>
00227 ThermalConductivity<EvalT,Traits>::getValidThermalCondParameters() const
00228 {
00229   Teuchos::RCP<Teuchos::ParameterList> validPL =
00230        rcp(new Teuchos::ParameterList("Valid Thermal Conductivity Params"));;
00231 
00232   validPL->set<std::string>("Thermal Conductivity Type", "Constant", 
00233                "Constant thermal conductivity across the entire domain");
00234   validPL->set<double>("Value", 1.0, "Constant thermal conductivity value");
00235 
00236 // Truncated KL parameters
00237 
00238   validPL->set<int>("Number of KL Terms", 2, "");
00239   validPL->set<double>("Mean", 0.2, "");
00240   validPL->set<double>("Standard Deviation", 0.1, "");
00241   validPL->set<std::string>("Domain Lower Bounds", "{0.0 0.0}", "");
00242   validPL->set<std::string>("Domain Upper Bounds", "{1.0 1.0}", "");
00243   validPL->set<std::string>("Correlation Lengths", "{1.0 1.0}", "");
00244   return validPL;
00245 }
00246 
00247 // **********************************************************************
00248 // **********************************************************************
00249 }
00250
PHAL_ThermalConductivity_Def.hpp
