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