00001
00002
00003
00004
00005
00006
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009
00010 #include "Intrepid_FunctionSpaceTools.hpp"
00011 #include "Intrepid_RealSpaceTools.hpp"
00012
00013 #include <typeinfo>
00014
00015 namespace LCM {
00016
00017
00018 template<typename EvalT, typename Traits>
00019 ThermoPoroPlasticityResidEnergy<EvalT, Traits>::
00020 ThermoPoroPlasticityResidEnergy(const Teuchos::ParameterList& p) :
00021 wBF (p.get<std::string> ("Weighted BF Name"),
00022 p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
00023 porePressure (p.get<std::string> ("QP Pore Pressure Name"),
00024 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00025 Temp (p.get<std::string> ("QP Temperature Name"),
00026 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00027 stabParameter (p.get<std::string> ("Material Property Name"),
00028 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00029 ThermalCond (p.get<std::string> ("Thermal Conductivity Name"),
00030 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00031 kcPermeability (p.get<std::string> ("Kozeny-Carman Permeability Name"),
00032 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00033 porosity (p.get<std::string> ("Porosity Name"),
00034 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00035 refTemp (p.get<std::string> ("Reference Temperature Name"),
00036 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00037 alphaSkeleton (p.get<std::string> ("Skeleton Thermal Expansion Name"),
00038 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00039 gammaMixture (p.get<std::string> ("Mixture Specific Heat Name"),
00040 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00041 gammaFluid (p.get<std::string> ("Pore-Fluid Specific Heat Name"),
00042 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00043 alphaMixture (p.get<std::string> ("Mixture Thermal Expansion Name"),
00044 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00045 biotCoefficient (p.get<std::string> ("Biot Coefficient Name"),
00046 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00047 biotModulus (p.get<std::string> ("Biot Modulus Name"),
00048 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00049 bulk (p.get<std::string> ("Bulk Modulus Name"),
00050 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00051 wGradBF (p.get<std::string> ("Weighted Gradient BF Name"),
00052 p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
00053 TGrad (p.get<std::string> ("Gradient QP Variable Name"),
00054 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00055 PGrad (p.get<std::string> ("Pore Pressure Gradient Name"),
00056 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00057 Source (p.get<std::string> ("Source Name"),
00058 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00059 coordVec (p.get<std::string> ("Coordinate Vector Name"),
00060 p.get<Teuchos::RCP<PHX::DataLayout> >("Coordinate Data Layout") ),
00061 cubature (p.get<Teuchos::RCP <Intrepid::Cubature<RealType> > >("Cubature")),
00062 cellType (p.get<Teuchos::RCP <shards::CellTopology> > ("Cell Type")),
00063 weights (p.get<std::string> ("Weights Name"),
00064 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00065 young_modulus_ (p.get<std::string> ("Elastic Modulus Name"),
00066 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00067 poissons_ratio_ (p.get<std::string> ("Poissons Ratio Name"),
00068 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00069 deltaTime (p.get<std::string>("Delta Time Name"),
00070 p.get<Teuchos::RCP<PHX::DataLayout> >("Workset Scalar Data Layout")),
00071 J (p.get<std::string> ("DetDefGrad Name"),
00072 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00073 defgrad (p.get<std::string> ("DefGrad Name"),
00074 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
00075 TResidual (p.get<std::string> ("Residual Name"),
00076 p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
00077 haveSource (p.get<bool>("Have Source")),
00078 haveConvection(false),
00079 haveAbsorption (p.get<bool>("Have Absorption")),
00080 haverhoCp(false)
00081 {
00082 if (p.isType<bool>("Disable Transient"))
00083 enableTransient = !p.get<bool>("Disable Transient");
00084 else enableTransient = true;
00085
00086 this->addDependentField(stabParameter);
00087 this->addDependentField(deltaTime);
00088 this->addDependentField(weights);
00089 this->addDependentField(coordVec);
00090 this->addDependentField(wBF);
00091 this->addDependentField(refTemp);
00092 this->addDependentField(alphaSkeleton);
00093 this->addDependentField(gammaMixture);
00094 this->addDependentField(gammaFluid);
00095 this->addDependentField(porePressure);
00096 this->addDependentField(ThermalCond);
00097 this->addDependentField(kcPermeability);
00098 this->addDependentField(porosity);
00099 this->addDependentField(biotCoefficient);
00100 this->addDependentField(biotModulus);
00101 this->addDependentField(young_modulus_);
00102 this->addDependentField(poissons_ratio_);
00103
00104 this->addDependentField(TGrad);
00105 this->addDependentField(PGrad);
00106 this->addDependentField(wGradBF);
00107 if (haveSource) this->addDependentField(Source);
00108 if (haveAbsorption) {
00109 Absorption = PHX::MDField<ScalarT,Cell,QuadPoint>(
00110 p.get<std::string>("Absorption Name"),
00111 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"));
00112 this->addDependentField(Absorption);
00113 }
00114
00115 this->addDependentField(J);
00116 this->addDependentField(defgrad);
00117 this->addDependentField(alphaMixture);
00118 this->addDependentField(bulk);
00119 this->addEvaluatedField(TResidual);
00120 this->addEvaluatedField(Temp);
00121
00122
00123
00124
00125
00126
00127
00128 Teuchos::RCP<PHX::DataLayout> vector_dl =
00129 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00130 std::vector<PHX::DataLayout::size_type> dims;
00131 vector_dl->dimensions(dims);
00132 numQPs = dims[1];
00133 numDims = dims[2];
00134
00135 Teuchos::RCP<PHX::DataLayout> node_dl =
00136 p.get< Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout");
00137 std::vector<PHX::DataLayout::size_type> ndims;
00138 node_dl->dimensions(ndims);
00139 worksetSize = dims[0];
00140 numNodes = dims[1];
00141
00142
00143
00144 porePressureName = p.get<std::string>("QP Pore Pressure Name")+"_old";
00145 JName =p.get<std::string>("DetDefGrad Name")+"_old";
00146 tempName =p.get<std::string>("QP Temperature Name")+"_old";
00147
00148
00149
00150
00151
00152
00153
00154
00155 C.resize(worksetSize, numQPs, numDims, numDims);
00156 Cinv.resize(worksetSize, numQPs, numDims, numDims);
00157 F_inv.resize(worksetSize, numQPs, numDims, numDims);
00158 F_invT.resize(worksetSize, numQPs, numDims, numDims);
00159 JF_invT.resize(worksetSize, numQPs, numDims, numDims);
00160 KJF_invT.resize(worksetSize, numQPs, numDims, numDims);
00161 Kref.resize(worksetSize, numQPs, numDims, numDims);
00162
00163
00164
00165
00166 flux.resize(dims[0], numQPs, numDims);
00167 fluxdt.resize(dims[0], numQPs, numDims);
00168 pterm.resize(dims[0], numQPs);
00169 tterm.resize(dims[0], numQPs);
00170
00171 tpterm.resize(dims[0], numNodes, numQPs);
00172
00173 if (haveAbsorption) aterm.resize(dims[0], numQPs);
00174
00175 convectionVels = Teuchos::getArrayFromStringParameter<double> (p,
00176 "Convection Velocity", numDims, false);
00177 if (p.isType<std::string>("Convection Velocity")) {
00178 convectionVels = Teuchos::getArrayFromStringParameter<double> (p,
00179 "Convection Velocity", numDims, false);
00180 }
00181 if (convectionVels.size()>0) {
00182 haveConvection = true;
00183 if (p.isType<bool>("Have Rho Cp"))
00184 haverhoCp = p.get<bool>("Have Rho Cp");
00185 if (haverhoCp) {
00186 PHX::MDField<ScalarT,Cell,QuadPoint> tmp(p.get<std::string>("Rho Cp Name"),
00187 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"));
00188 rhoCp = tmp;
00189 this->addDependentField(rhoCp);
00190 }
00191 }
00192
00193 this->setName("ThermoPoroPlasticityResidEnergy"+PHX::TypeString<EvalT>::value);
00194
00195 }
00196
00197
00198 template<typename EvalT, typename Traits>
00199 void ThermoPoroPlasticityResidEnergy<EvalT, Traits>::
00200 postRegistrationSetup(typename Traits::SetupData d,
00201 PHX::FieldManager<Traits>& fm)
00202 {
00203 this->utils.setFieldData(stabParameter,fm);
00204 this->utils.setFieldData(deltaTime,fm);
00205 this->utils.setFieldData(weights,fm);
00206 this->utils.setFieldData(coordVec,fm);
00207 this->utils.setFieldData(wBF,fm);
00208 this->utils.setFieldData(porePressure,fm);
00209 this->utils.setFieldData(ThermalCond,fm);
00210 this->utils.setFieldData(kcPermeability,fm);
00211 this->utils.setFieldData(porosity,fm);
00212 this->utils.setFieldData(alphaMixture,fm);
00213 this->utils.setFieldData(biotCoefficient,fm);
00214 this->utils.setFieldData(biotModulus,fm);
00215 this->utils.setFieldData(bulk,fm);
00216 this->utils.setFieldData(TGrad,fm);
00217 this->utils.setFieldData(PGrad,fm);
00218 this->utils.setFieldData(wGradBF,fm);
00219 if (haveSource) this->utils.setFieldData(Source,fm);
00220 this->utils.setFieldData(Temp,fm);
00221 if (haveAbsorption) this->utils.setFieldData(Absorption,fm);
00222 if (haveConvection && haverhoCp) this->utils.setFieldData(rhoCp,fm);
00223 this->utils.setFieldData(J,fm);
00224 this->utils.setFieldData(defgrad,fm);
00225 this->utils.setFieldData(refTemp,fm);
00226 this->utils.setFieldData(alphaSkeleton,fm);
00227 this->utils.setFieldData(gammaMixture,fm);
00228 this->utils.setFieldData(gammaFluid,fm);
00229 this->utils.setFieldData(young_modulus_,fm);
00230 this->utils.setFieldData(poissons_ratio_,fm);
00231 this->utils.setFieldData(TResidual,fm);
00232 }
00233
00234
00235 template<typename EvalT, typename Traits>
00236 void ThermoPoroPlasticityResidEnergy<EvalT, Traits>::
00237 evaluateFields(typename Traits::EvalData workset)
00238 {
00239 typedef Intrepid::FunctionSpaceTools FST;
00240 typedef Intrepid::RealSpaceTools<ScalarT> RST;
00241
00242 Albany::MDArray porePressureold = (*workset.stateArrayPtr)[porePressureName];
00243 Albany::MDArray Jold = (*workset.stateArrayPtr)[JName];
00244 Albany::MDArray Tempold = (*workset.stateArrayPtr)[tempName];
00245
00246 ScalarT dTemperature(0.0);
00247 ScalarT dporePressure(0.0);
00248 ScalarT dJ(0.0);
00249
00250
00251
00252 ScalarT dt = deltaTime(0);
00253
00254 RST::inverse(F_inv, defgrad);
00255 RST::transpose(F_invT, F_inv);
00256 FST::scalarMultiplyDataData<ScalarT>(JF_invT, J, F_invT);
00257 FST::scalarMultiplyDataData<ScalarT>(KJF_invT, ThermalCond, JF_invT);
00258 FST::tensorMultiplyDataData<ScalarT>(Kref, F_inv, KJF_invT);
00259
00260 FST::tensorMultiplyDataData<ScalarT> (flux, Kref, TGrad);
00261
00262 for (std::size_t cell=0; cell < workset.numCells; ++cell){
00263 for (std::size_t qp=0; qp < numQPs; ++qp) {
00264 for (std::size_t dim=0; dim <numDims; ++dim){
00265 fluxdt(cell, qp, dim) = flux(cell,qp,dim)*dt;
00266 }
00267 }
00268 }
00269
00270 FST::integrate<ScalarT>(TResidual, fluxdt, wGradBF, Intrepid::COMP_CPP, false);
00271
00272
00273 FST::scalarMultiplyDataData<ScalarT>(KJF_invT, kcPermeability, JF_invT);
00274 FST::tensorMultiplyDataData<ScalarT>(Kref, F_inv, KJF_invT);
00275 FST::tensorMultiplyDataData<ScalarT> (flux, Kref, PGrad);
00276 FST::tensorMultiplyDataData<ScalarT> (fluxdt, F_invT, flux);
00277
00278
00279 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00280
00281 for (std::size_t node=0; node < numNodes; ++node) {
00282 for (std::size_t qp=0; qp < numQPs; ++qp) {
00283 for (std::size_t dim=0; dim <numDims; ++dim){
00284 TResidual(cell,node) -= dt*gammaFluid(cell,qp)*
00285 porosity(cell,qp)*
00286 fluxdt(cell,qp,dim)*
00287 TGrad(cell,qp, dim)*
00288 wBF(cell, node, qp) ;
00289 }
00290 }
00291 }
00292 }
00293
00294
00295
00296
00297 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00298
00299 for (std::size_t node=0; node < numNodes; ++node) {
00300
00301 for (std::size_t qp=0; qp < numQPs; ++qp) {
00302
00303
00304 dJ = std::log(J(cell,qp)/Jold(cell,qp));
00305 dTemperature = Temp(cell,qp)-Tempold(cell,qp);
00306 dporePressure = porePressure(cell,qp)-porePressureold(cell, qp);
00307
00308
00309 TResidual(cell,node) += 3.0*alphaSkeleton(cell,qp)
00310 *bulk(cell,qp)*Temp(cell,qp)
00311 *dJ*wBF(cell, node, qp) ;
00312
00313
00314 TResidual(cell,node) -= 3.0*dporePressure
00315 *alphaMixture(cell,qp)
00316 *Temp(cell,qp)
00317 *wBF(cell, node, qp);
00318
00319 TResidual(cell,node) += dTemperature
00320 *gammaMixture(cell, qp)
00321 *wBF(cell, node, qp);
00322
00323 }
00324 }
00325 }
00326
00327
00328
00329
00330
00331
00332
00333 for (std::size_t cell=0; cell < workset.numCells; ++cell){
00334
00335 porePbar = 0.0;
00336 Tempbar = 0.0;
00337 vol = 0.0;
00338 for (std::size_t qp=0; qp < numQPs; ++qp) {
00339
00340 porePbar += weights(cell,qp)*(
00341 (porePressure(cell,qp)-porePressureold(cell, qp) ));
00342 Tempbar += weights(cell,qp)*(Temp(cell,qp)-Tempold(cell,qp));
00343
00344
00345 vol += weights(cell,qp);
00346 }
00347 porePbar /= vol;
00348 Tempbar /= vol;
00349
00350 for (std::size_t qp=0; qp < numQPs; ++qp) {
00351 pterm(cell,qp) = porePbar;
00352 tterm(cell,qp) = Tempbar;
00353
00354 }
00355 }
00356
00357 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00358 for (std::size_t node=0; node < numNodes; ++node) {
00359 for (std::size_t qp=0; qp < numQPs; ++qp) {
00360
00361 shearModulus = young_modulus_(cell,qp)*0.5/(1.0+ poissons_ratio_(cell,qp));
00362 bulkModulus = young_modulus_(cell,qp)/3.0/(1.0- 2.0*poissons_ratio_(cell,qp));
00363
00364 dTemperature = Temp(cell,qp)-Tempold(cell,qp) -tterm(cell,qp);
00365 dporePressure = porePressure(cell,qp)-porePressureold(cell, qp)- pterm(cell,qp);
00366
00367 TResidual(cell,node) -= 3.0*dporePressure
00368 *Temp(cell,qp)
00369 *stabParameter(cell, qp)
00370 *alphaMixture(cell,qp)
00371 *wBF(cell, node, qp);
00372
00373 TResidual(cell,node) += (dTemperature)
00374 *stabParameter(cell, qp)
00375 *(gammaMixture(cell, qp) + Temp(cell,qp)
00376 *9.0*alphaSkeleton(cell,qp)*bulkModulus
00377 *alphaSkeleton(cell,qp)*bulkModulus
00378 /(bulkModulus + 4.0/3.0*shearModulus))
00379 *wBF(cell, node, qp);
00380
00381 }
00382 }
00383 }
00384 }
00385
00386 }
00387
00388
