Albany: a Trilinos-based PDE code: ConstitutiveModelParameters

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 "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 #include "Sacado_ParameterRegistration.hpp"
00010 #include "Albany_Utils.hpp"
00011 #include "Teuchos_Array.hpp"
00012 
00013 namespace LCM
00014 {
00015 
00016 //------------------------------------------------------------------------------
00017 template<typename EvalT, typename Traits>
00018 ConstitutiveModelParameters<EvalT, Traits>::
00019 ConstitutiveModelParameters(Teuchos::ParameterList& p,
00020     const Teuchos::RCP<Albany::Layouts>& dl) :
00021     have_temperature_(false),
00022     dl_(dl)
00023 {
00024   // get number of integration points and spatial dimensions
00025   std::vector<PHX::DataLayout::size_type> dims;
00026   dl_->qp_vector->dimensions(dims);
00027   num_pts_ = dims[1];
00028   num_dims_ = dims[2];
00029 
00030   // get the Parameter Library
00031   Teuchos::RCP<ParamLib> paramLib =
00032       p.get<Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00033 
00034   // get the material parameter list
00035   Teuchos::ParameterList* mat_params =
00036       p.get<Teuchos::ParameterList*>("Material Parameters");
00037 
00038   // Check for optional field: temperature
00039   if (p.isType<std::string>("Temperature Name")) {
00040     have_temperature_ = true;
00041     PHX::MDField<ScalarT, Cell, QuadPoint>
00042     tmp(p.get<std::string>("Temperature Name"), dl_->qp_scalar);
00043     temperature_ = tmp;
00044     this->addDependentField(temperature_);
00045   }
00046 
00047   // step through the possible parameters, registering as necessary
00048   //
00049   // elastic modulus
00050   std::string e_mod("Elastic Modulus");
00051   if (mat_params->isSublist(e_mod)) {
00052     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(e_mod, dl_->qp_scalar);
00053     elastic_mod_ = tmp;
00054     field_map_.insert(std::make_pair(e_mod, elastic_mod_));
00055     parseParameters(e_mod, p, paramLib);
00056   }
00057   // Poisson's ratio
00058   std::string pr("Poissons Ratio");
00059   if (mat_params->isSublist(pr)) {
00060     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(pr, dl_->qp_scalar);
00061     poissons_ratio_ = tmp;
00062     field_map_.insert(std::make_pair(pr, poissons_ratio_));
00063     parseParameters(pr, p, paramLib);
00064   }
00065   // bulk modulus
00066   std::string b_mod("Bulk Modulus");
00067   if (mat_params->isSublist(b_mod)) {
00068     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(b_mod, dl_->qp_scalar);
00069     bulk_mod_ = tmp;
00070     field_map_.insert(std::make_pair(b_mod, bulk_mod_));
00071     parseParameters(b_mod, p, paramLib);
00072   }
00073   // shear modulus
00074   std::string s_mod("Shear Modulus");
00075   if (mat_params->isSublist(s_mod)) {
00076     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(s_mod, dl_->qp_scalar);
00077     shear_mod_ = tmp;
00078     field_map_.insert(std::make_pair(s_mod, shear_mod_));
00079     parseParameters(s_mod, p, paramLib);
00080   }
00081   // yield strength
00082   std::string yield("Yield Strength");
00083   if (mat_params->isSublist(yield)) {
00084     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(yield, dl_->qp_scalar);
00085     yield_strength_ = tmp;
00086     field_map_.insert(std::make_pair(yield, yield_strength_));
00087     parseParameters(yield, p, paramLib);
00088   }
00089   // hardening modulus
00090   std::string h_mod("Hardening Modulus");
00091   if (mat_params->isSublist(h_mod)) {
00092     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(h_mod, dl_->qp_scalar);
00093     hardening_mod_ = tmp;
00094     field_map_.insert(std::make_pair(h_mod, hardening_mod_));
00095     parseParameters(h_mod, p, paramLib);
00096   }
00097   // recovery modulus
00098   std::string r_mod("Recovery Modulus");
00099   if (mat_params->isSublist(r_mod)) {
00100     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(r_mod, dl_->qp_scalar);
00101     recovery_mod_ = tmp;
00102     field_map_.insert(std::make_pair(r_mod, recovery_mod_));
00103     parseParameters(r_mod, p, paramLib);
00104   }
00105   // concentration equilibrium parameter
00106   std::string c_eq("Concentration Equilibrium Parameter");
00107   if (mat_params->isSublist(c_eq)) {
00108     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(c_eq, dl_->qp_scalar);
00109     conc_eq_param_ = tmp;
00110     field_map_.insert(std::make_pair(c_eq, conc_eq_param_));
00111     parseParameters(c_eq, p, paramLib);
00112   }
00113   // diffusion coefficient
00114   std::string d_coeff("Diffusion Coefficient");
00115   if (mat_params->isSublist(d_coeff)) {
00116     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(d_coeff, dl_->qp_scalar);
00117     diff_coeff_ = tmp;
00118     field_map_.insert(std::make_pair(d_coeff, diff_coeff_));
00119     parseParameters(d_coeff, p, paramLib);
00120   }
00121   // thermal conductivity
00122   std::string th_cond("Thermal Conductivity");
00123   if (mat_params->isSublist(th_cond)) {
00124     PHX::MDField<ScalarT, Cell, QuadPoint> tmp(th_cond, dl_->qp_scalar);
00125     thermal_cond_ = tmp;
00126     field_map_.insert(std::make_pair(th_cond, thermal_cond_));
00127     parseParameters(th_cond, p, paramLib);
00128   }
00129 
00130   // register evaluated fields
00131   typename
00132   std::map<std::string, PHX::MDField<ScalarT, Cell, QuadPoint> >::iterator it;
00133   for (it = field_map_.begin();
00134       it != field_map_.end();
00135       ++it) {
00136     this->addEvaluatedField(it->second);
00137   }
00138   this->setName(
00139       "Constitutive Model Parameters" + PHX::TypeString<EvalT>::value);
00140 }
00141 
00142 //------------------------------------------------------------------------------
00143 template<typename EvalT, typename Traits>
00144 void ConstitutiveModelParameters<EvalT, Traits>::
00145 postRegistrationSetup(typename Traits::SetupData d,
00146     PHX::FieldManager<Traits>& fm)
00147 {
00148   typename std::map<std::string, PHX::MDField<ScalarT, Cell, QuadPoint> >::iterator it;
00149   for (it = field_map_.begin();
00150       it != field_map_.end();
00151       ++it) {
00152     this->utils.setFieldData(it->second, fm);
00153     if (!is_constant_map_[it->first]) {
00154       this->utils.setFieldData(coord_vec_, fm);
00155     }
00156   }
00157 
00158   if (have_temperature_) this->utils.setFieldData(temperature_, fm);
00159 }
00160 
00161 //------------------------------------------------------------------------------
00162 template<typename EvalT, typename Traits>
00163 void ConstitutiveModelParameters<EvalT, Traits>::
00164 evaluateFields(typename Traits::EvalData workset)
00165 {
00166   typename std::map<std::string, PHX::MDField<ScalarT, Cell, QuadPoint> >::iterator it;
00167   for (it = field_map_.begin();
00168       it != field_map_.end();
00169       ++it) {
00170     ScalarT constant_value = constant_value_map_[it->first];
00171     if (is_constant_map_[it->first]) {
00172       for (std::size_t cell(0); cell < workset.numCells; ++cell) {
00173         for (std::size_t pt(0); pt < num_pts_; ++pt) {
00174           it->second(cell, pt) = constant_value;
00175         }
00176       }
00177     } else {
00178       for (std::size_t cell(0); cell < workset.numCells; ++cell) {
00179         for (std::size_t pt(0); pt < num_pts_; ++pt) {
00180           Teuchos::Array<MeshScalarT> point(num_dims_);
00181           for (std::size_t i(0); i < num_dims_; ++i)
00182             point[i] = Sacado::ScalarValue<MeshScalarT>::eval(
00183                 coord_vec_(cell, pt, i));
00184           it->second(cell, pt) =
00185               exp_rf_kl_map_[it->first]->evaluate(point, rv_map_[it->first]);
00186         }
00187       }
00188     }
00189     // FIXME deal with Arrhenius temperature dependence, too
00190     if (have_temperature_) {
00191       RealType dPdT = dparam_dtemp_map_[it->first];
00192       RealType ref_temp = ref_temp_map_[it->first];
00193       for (std::size_t cell(0); cell < workset.numCells; ++cell) {
00194         for (std::size_t pt(0); pt < num_pts_; ++pt) {
00195           it->second(cell, pt) += dPdT * (temperature_(cell, pt) - ref_temp);
00196         }
00197       }
00198     }
00199   }
00200 }
00201 
00202 //------------------------------------------------------------------------------
00203 template<typename EvalT, typename Traits>
00204 typename ConstitutiveModelParameters<EvalT, Traits>::ScalarT&
00205 ConstitutiveModelParameters<EvalT, Traits>::getValue(const std::string &n)
00206 {
00207   typename std::map<std::string, ScalarT>::iterator it;
00208   for (it = constant_value_map_.begin();
00209       it != constant_value_map_.end();
00210       ++it) {
00211     if (n == it->first) {
00212       return constant_value_map_[it->first];
00213     }
00214   }
00215   typename std::map<std::string, Teuchos::Array<ScalarT> >::iterator it2;
00216   for (int i(0); i < rv_map_[it2->first].size(); ++i) {
00217     if (n == Albany::strint(n + " KL Random Variable", i))
00218       return rv_map_[it2->first][i];
00219   }
00220 }
00221 //------------------------------------------------------------------------------
00222 template<typename EvalT, typename Traits>
00223 void ConstitutiveModelParameters<EvalT, Traits>::
00224 parseParameters(const std::string &n,
00225     Teuchos::ParameterList &p,
00226     Teuchos::RCP<ParamLib> paramLib)
00227 {
00228   Teuchos::ParameterList pl =
00229     p.get<Teuchos::ParameterList*>("Material Parameters")->sublist(n);
00230   std::string type_name(n + " Type");
00231   std::string type = pl.get(type_name, "Constant");
00232   if (type == "Constant") {
00233     is_constant_map_.insert(std::make_pair(n, true));
00234     constant_value_map_.insert(std::make_pair(n, pl.get("Value", 1.0)));
00235     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(n, this, paramLib);
00236     if (have_temperature_) {
00237       if (pl.get<std::string>("Temperature Dependence Type", "Linear")
00238           == "Linear") {
00239         dparam_dtemp_map_.insert
00240         (std::make_pair(n,
00241           pl.get<RealType>("Linear Temperature Coefficient", 0.0)));
00242         ref_temp_map_.insert
00243         (std::make_pair(n, pl.get<RealType>("Reference Temperature", -1)));
00244       } else if (pl.get<std::string>("Temperature Dependence Type", "Linear")
00245           == "Arrhenius") {
00246         ideal_map_.insert(
00247             std::make_pair(n, pl.get<RealType>("Ideal Gas Constant", 1.0)));
00248         pre_exp_map_.insert(
00249             std::make_pair(n, pl.get<RealType>("Pre Exponential", 0.0)));
00250         exp_param_map_.insert(
00251             std::make_pair(n, pl.get<RealType>("Exponential Parameter", 0.0)));
00252       }
00253     }
00254   } else if (type == "Truncated KL Expansion") {
00255     is_constant_map_.insert(std::make_pair(n, false));
00256     PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim>
00257       fx(p.get<std::string>("QP Coordinate Vector Name"), dl_->qp_vector);
00258     coord_vec_ = fx;
00259     this->addDependentField(coord_vec_);
00260 
00261     exp_rf_kl_map_.
00262         insert(
00263         std::make_pair(n,
00264             Teuchos::rcp(
00265                 new Stokhos::KL::ExponentialRandomField<MeshScalarT>(pl))));
00266     int num_KL = exp_rf_kl_map_[n]->stochasticDimension();
00267 
00268     // Add KL random variables as Sacado-ized parameters
00269     rv_map_.insert(std::make_pair(n, Teuchos::Array<ScalarT>(num_KL)));
00270     for (int i(0); i < num_KL; ++i) {
00271       std::string ss = Albany::strint(n + " KL Random Variable", i);
00272       new Sacado::ParameterRegistration<EvalT, SPL_Traits>(ss, this, paramLib);
00273       rv_map_[n][i] = pl.get(ss, 0.0);
00274     }
00275   }
00276 }
00277 //------------------------------------------------------------------------------
00278 }
00279
ConstitutiveModelParameters_Def.hpp
