Albany: a Trilinos-based PDE code: TLPoroPlasticityResidMass

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 
00010 #include "Intrepid_FunctionSpaceTools.hpp"
00011 #include "Intrepid_RealSpaceTools.hpp"
00012 
00013 #include <typeinfo>
00014 namespace LCM {
00015 
00016   //**********************************************************************
00017   template<typename EvalT, typename Traits>
00018   TLPoroPlasticityResidMass<EvalT, Traits>::
00019   TLPoroPlasticityResidMass(Teuchos::ParameterList& p) :
00020     wBF   (p.get<std::string>                   ("Weighted BF Name"),
00021          p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
00022     porePressure (p.get<std::string>                   ("QP Pore Pressure Name"),
00023          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00024     elementLength (p.get<std::string>         ("Element Length Name"),
00025          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00026     Tdot   (p.get<std::string>                   ("QP Time Derivative Variable Name"),
00027      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00028     ThermalCond (p.get<std::string>                   ("Thermal Conductivity Name"),
00029      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00030     kcPermeability (p.get<std::string>            ("Kozeny-Carman Permeability Name"),
00031          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00032     porosity (p.get<std::string>                   ("Porosity Name"),
00033          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00034     biotCoefficient (p.get<std::string>           ("Biot Coefficient Name"),
00035          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00036     biotModulus (p.get<std::string>                   ("Biot Modulus Name"),
00037      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00038     wGradBF  (p.get<std::string>                   ("Weighted Gradient BF Name"),
00039      p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
00040     TGrad       (p.get<std::string>                   ("Gradient QP Variable Name"),
00041      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00042     Source      (p.get<std::string>                   ("Source Name"),
00043      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00044     coordVec  (p.get<std::string>                   ("Coordinate Vector Name"),
00045          p.get<Teuchos::RCP<PHX::DataLayout> >("Coordinate Data Layout") ),
00046     cubature   (p.get<Teuchos::RCP <Intrepid::Cubature<RealType> > >("Cubature")),
00047     cellType    (p.get<Teuchos::RCP <shards::CellTopology> > ("Cell Type")),
00048     weights     (p.get<std::string>                   ("Weights Name"),
00049          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00050     deltaTime (p.get<std::string>("Delta Time Name"),
00051          p.get<Teuchos::RCP<PHX::DataLayout> >("Workset Scalar Data Layout")),
00052     TResidual  (p.get<std::string>                   ("Residual Name"),
00053      p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
00054     haveSource  (p.get<bool>("Have Source")),
00055     haveConvection(false),
00056     haveAbsorption(p.get<bool>("Have Absorption")),
00057     haverhoCp(false),
00058     haveMechanics(p.get<bool>("Have Mechanics", false)),
00059     stab_param_(p.get<RealType>("Stabilization Parameter"))
00060   {
00061     //  if (p.isType<bool>("Disable Transient"))
00062     //    enableTransient = !p.get<bool>("Disable Transient");
00063     //  else enableTransient = true;
00064 
00065     enableTransient = false;
00066 
00067       // this->addDependentField(stabParameter);
00068     this->addDependentField(elementLength);
00069     this->addDependentField(deltaTime);
00070     this->addDependentField(weights);
00071     this->addDependentField(coordVec);
00072     this->addDependentField(wBF);
00073     this->addDependentField(porePressure);
00074     this->addDependentField(ThermalCond);
00075     this->addDependentField(kcPermeability);
00076     this->addDependentField(porosity);
00077     this->addDependentField(biotCoefficient);
00078     this->addDependentField(biotModulus);
00079     if (enableTransient) this->addDependentField(Tdot);
00080     this->addDependentField(TGrad);
00081     this->addDependentField(wGradBF);
00082     if (haveSource) this->addDependentField(Source);
00083     if (haveAbsorption) {
00084       Absorption = PHX::MDField<ScalarT,Cell,QuadPoint>
00085         (p.get<std::string>("Absorption Name"),
00086          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"));
00087       this->addDependentField(Absorption);
00088     }
00089 
00090     if (p.isType<std::string>("DefGrad Name")) {
00091       Teuchos::RCP<PHX::DataLayout> tensor_dl =
00092         p.get< Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout");
00093       Teuchos::RCP<PHX::DataLayout> scalar_dl =
00094         p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00095 
00096       haveMechanics = true;
00097 
00098       PHX::MDField<ScalarT,Cell,QuadPoint,Dim, Dim>
00099         tf(p.get<std::string>("DefGrad Name"), tensor_dl);
00100       defgrad = tf;
00101       this->addDependentField(defgrad);
00102 
00103       PHX::MDField<ScalarT,Cell,QuadPoint>
00104         tj(p.get<std::string>("DetDefGrad Name"), scalar_dl);
00105       J = tj;
00106       this->addDependentField(J);
00107     }
00108 
00109     this->addEvaluatedField(TResidual);
00110 
00111 
00112     Teuchos::RCP<PHX::DataLayout> vector_dl =
00113       p.get< Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout");
00114     std::vector<PHX::DataLayout::size_type> dims;
00115     vector_dl->dimensions(dims);
00116 
00117     // Get data from previous converged time step
00118     porePressureName = p.get<std::string>("QP Pore Pressure Name")+"_old";
00119     if (haveMechanics) JName = p.get<std::string>("DetDefGrad Name")+"_old";
00120 
00121     worksetSize = dims[0];
00122     numNodes = dims[1];
00123     numQPs  = dims[2];
00124     numDims = dims[3];
00125 
00126     if (haveMechanics) {
00127       // Works space FCs
00128       C.resize(worksetSize, numQPs, numDims, numDims);
00129       Cinv.resize(worksetSize, numQPs, numDims, numDims);
00130       F_inv.resize(worksetSize, numQPs, numDims, numDims);
00131       F_invT.resize(worksetSize, numQPs, numDims, numDims);
00132       JF_invT.resize(worksetSize, numQPs, numDims, numDims);
00133       KJF_invT.resize(worksetSize, numQPs, numDims, numDims);
00134       Kref.resize(worksetSize, numQPs, numDims, numDims);
00135     }
00136 
00137     // Allocate workspace
00138     flux.resize(dims[0], numQPs, numDims);
00139     fluxdt.resize(dims[0], numQPs, numDims);
00140     pterm.resize(dims[0], numQPs);
00141     tpterm.resize(dims[0], numNodes, numQPs);
00142 
00143     if (haveAbsorption)  aterm.resize(dims[0], numQPs);
00144 
00145     convectionVels = Teuchos::getArrayFromStringParameter<double> 
00146       (p,"Convection Velocity",numDims,false);
00147     if (p.isType<std::string>("Convection Velocity")) {
00148       convectionVels = Teuchos::getArrayFromStringParameter<double> 
00149         (p,"Convection Velocity",numDims,false);
00150     }
00151     if (convectionVels.size()>0) {
00152       haveConvection = true;
00153       if (p.isType<bool>("Have Rho Cp"))
00154         haverhoCp = p.get<bool>("Have Rho Cp");
00155       if (haverhoCp) {
00156         PHX::MDField<ScalarT,Cell,QuadPoint> 
00157           tmp(p.get<std::string>("Rho Cp Name"),
00158               p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"));
00159         rhoCp = tmp;
00160         this->addDependentField(rhoCp);
00161       }
00162     }
00163 
00164     this->setName("TLPoroPlasticityResidMass"+PHX::TypeString<EvalT>::value);
00165 
00166   }
00167 
00168   //**********************************************************************
00169   template<typename EvalT, typename Traits>
00170   void TLPoroPlasticityResidMass<EvalT, Traits>::
00171   postRegistrationSetup(typename Traits::SetupData d,
00172                         PHX::FieldManager<Traits>& fm)
00173   {
00174     //this->utils.setFieldData(stabParameter,fm);
00175     this->utils.setFieldData(elementLength,fm);
00176     this->utils.setFieldData(deltaTime,fm);
00177     this->utils.setFieldData(weights,fm);
00178     this->utils.setFieldData(coordVec,fm);
00179     this->utils.setFieldData(wBF,fm);
00180     this->utils.setFieldData(porePressure,fm);
00181     this->utils.setFieldData(ThermalCond,fm);
00182     this->utils.setFieldData(kcPermeability,fm);
00183     this->utils.setFieldData(porosity,fm);
00184     this->utils.setFieldData(biotCoefficient,fm);
00185     this->utils.setFieldData(biotModulus,fm);
00186     this->utils.setFieldData(TGrad,fm);
00187     this->utils.setFieldData(wGradBF,fm);
00188     if (haveSource)  this->utils.setFieldData(Source,fm);
00189     if (enableTransient) this->utils.setFieldData(Tdot,fm);
00190     if (haveAbsorption)  this->utils.setFieldData(Absorption,fm);
00191     if (haveConvection && haverhoCp)  this->utils.setFieldData(rhoCp,fm);
00192     if (haveMechanics) {
00193       this->utils.setFieldData(J,fm);
00194       this->utils.setFieldData(defgrad,fm);
00195     }
00196     this->utils.setFieldData(TResidual,fm);
00197   }
00198 
00199   //**********************************************************************
00200   template<typename EvalT, typename Traits>
00201   void TLPoroPlasticityResidMass<EvalT, Traits>::
00202   evaluateFields(typename Traits::EvalData workset)
00203   {
00204     bool print = false;
00205     //if (typeid(ScalarT) == typeid(RealType)) print = true;
00206 
00207     typedef Intrepid::FunctionSpaceTools FST;
00208     typedef Intrepid::RealSpaceTools<ScalarT> RST;
00209 
00210     // Use previous time step for Backward Euler Integration
00211     Albany::MDArray porePressureold
00212       = (*workset.stateArrayPtr)[porePressureName];
00213 
00214     Albany::MDArray Jold;
00215     if (haveMechanics) {
00216       Jold = (*workset.stateArrayPtr)[JName];
00217     } 
00218 
00219     // Pore-fluid diffusion coupling.
00220     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00221       for (std::size_t node=0; node < numNodes; ++node) {
00222         TResidual(cell,node)=0.0;
00223         for (std::size_t qp=0; qp < numQPs; ++qp) {
00224 
00225           // Volumetric Constraint Term
00226           if (haveMechanics)  {
00227             TResidual(cell,node) -= biotCoefficient(cell, qp)
00228               * (std::log(J(cell,qp)/Jold(cell,qp)))
00229               * wBF(cell, node, qp) ;
00230           }
00231 
00232           // Pore-fluid Resistance Term
00233           TResidual(cell,node) -=  ( porePressure(cell,qp)-porePressureold(cell, qp) )
00234             / biotModulus(cell, qp)*wBF(cell, node, qp);
00235 
00236         }
00237       }
00238     }
00239     // Pore-Fluid Diffusion Term
00240 
00241     ScalarT dt = deltaTime(0);
00242 
00243     if (haveMechanics) {
00244       RST::inverse(F_inv, defgrad);
00245       RST::transpose(F_invT, F_inv);
00246       FST::scalarMultiplyDataData<ScalarT>(JF_invT, J, F_invT);
00247       FST::scalarMultiplyDataData<ScalarT>(KJF_invT, kcPermeability, JF_invT);
00248       FST::tensorMultiplyDataData<ScalarT>(Kref, F_inv, KJF_invT);
00249       FST::tensorMultiplyDataData<ScalarT> (flux, Kref, TGrad); // flux_i = k I_ij p_j
00250     } else {
00251       FST::scalarMultiplyDataData<ScalarT> (flux, kcPermeability, TGrad); // flux_i = kc p_i
00252     }
00253 
00254     for (std::size_t cell=0; cell < workset.numCells; ++cell){
00255       for (std::size_t qp=0; qp < numQPs; ++qp) {
00256         for (std::size_t dim=0; dim <numDims; ++dim){
00257           fluxdt(cell, qp, dim) = -flux(cell,qp,dim)*dt;
00258         }
00259       }
00260     }
00261     FST::integrate<ScalarT>(TResidual, fluxdt,
00262                             wGradBF, Intrepid::COMP_CPP, true); // "true" sums into
00263 
00264     //---------------------------------------------------------------------------//
00265     // Stabilization Term
00266 
00267     for (std::size_t cell=0; cell < workset.numCells; ++cell){
00268 
00269       porePbar = 0.0;
00270       vol = 0.0;
00271       for (std::size_t qp=0; qp < numQPs; ++qp) {
00272         porePbar += weights(cell,qp)
00273           * (porePressure(cell,qp)-porePressureold(cell, qp) );
00274         vol  += weights(cell,qp);
00275       }
00276       porePbar /= vol;
00277       for (std::size_t qp=0; qp < numQPs; ++qp) {
00278         pterm(cell,qp) = porePbar;
00279       }
00280 
00281       for (std::size_t node=0; node < numNodes; ++node) {
00282         trialPbar = 0.0;
00283         for (std::size_t qp=0; qp < numQPs; ++qp) {
00284           trialPbar += wBF(cell,node,qp);
00285         }
00286         trialPbar /= vol;
00287         for (std::size_t qp=0; qp < numQPs; ++qp) {
00288           tpterm(cell,node,qp) = trialPbar;
00289         }
00290 
00291       }
00292 
00293     }
00294     ScalarT temp(0);
00295 
00296     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00297 
00298       for (std::size_t node=0; node < numNodes; ++node) {
00299         for (std::size_t qp=0; qp < numQPs; ++qp) {
00300 
00301           temp = 3.0 - 12.0*kcPermeability(cell,qp)*dt
00302             /(elementLength(cell,qp)*elementLength(cell,qp));
00303 
00304           //if ((temp > 0) & stabParameter(cell,qp) > 0) {
00305           if ((temp > 0) & stab_param_ > 0) {
00306 
00307             TResidual(cell,node) -= 
00308               ( porePressure(cell,qp)-porePressureold(cell, qp) )
00309               //* stabParameter(cell, qp)
00310               * stab_param_
00311               * std::abs(temp) // should be 1 but use 0.5 for safety
00312               * (0.5 + 0.5*std::tanh( (temp-1)/kcPermeability(cell,qp)  ))
00313               / biotModulus(cell, qp)
00314               * ( wBF(cell, node, qp)
00315                 // -tpterm(cell,node,qp)
00316                 );
00317             TResidual(cell,node) += pterm(cell,qp)
00318               //* stabParameter(cell, qp)
00319               * stab_param_
00320               * std::abs(temp) // should be 1 but use 0.5 for safety
00321               * (0.5 + 0.5*std::tanh( (temp-1)/kcPermeability(cell,qp)  ))
00322               / biotModulus(cell, qp)
00323               * ( wBF(cell, node, qp) );
00324           }
00325         }
00326       }
00327     }
00328   }
00329   //**********************************************************************
00330 }
TLPoroPlasticityResidMass_Def.hpp
