Albany: a Trilinos-based PDE code: ThermoPoroPlasticityResidMass

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 "Intrepid_FunctionSpaceTools.hpp"
00010 #include "Intrepid_RealSpaceTools.hpp"
00011 #include <Intrepid_MiniTensor.h>
00012 
00013 #include <typeinfo>
00014 
00015 namespace LCM {
00016 
00017   //**********************************************************************
00018   template<typename EvalT, typename Traits>
00019   ThermoPoroPlasticityResidMass<EvalT, Traits>::
00020   ThermoPoroPlasticityResidMass(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   densityPoreFluid       (p.get<std::string>      ("Pore-Fluid Density Name"),
00026        p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00027     Temp        (p.get<std::string>                   ("QP Temperature Name"),
00028      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00029     RefTemp        (p.get<std::string>                   ("Reference Temperature Name"),
00030      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00031   young_modulus_ (p.get<std::string>                   ("Elastic Modulus Name"),
00032          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00033   poissons_ratio_       (p.get<std::string>      ("Poissons Ratio Name"),
00034            p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00035   stabParameter        (p.get<std::string>           ("Material Property Name"),
00036      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00037     ThermalCond (p.get<std::string>                   ("Thermal Conductivity Name"),
00038      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00039     kcPermeability (p.get<std::string>            ("Kozeny-Carman Permeability Name"),
00040      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00041     porosity (p.get<std::string>                   ("Porosity 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     alphaPoreFluid       (p.get<std::string>      ("Pore-Fluid Thermal Expansion Name"),
00046        p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00047   alphaSkeleton       (p.get<std::string>      ("Skeleton Thermal Expansion Name"),
00048        p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00049     biotCoefficient (p.get<std::string>           ("Biot Coefficient Name"),
00050      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00051     biotModulus (p.get<std::string>                   ("Biot Modulus Name"),
00052      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00053     wGradBF     (p.get<std::string>                   ("Weighted Gradient BF Name"),
00054      p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
00055     TGrad       (p.get<std::string>                   ("Gradient QP Variable Name"),
00056      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00057   TempGrad       (p.get<std::string>                   ("Temperature Gradient Name"),
00058      p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00059   weights       (p.get<std::string>                   ("Weights Name"),
00060     p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00061   deltaTime (p.get<std::string>("Delta Time Name"),
00062     p.get<Teuchos::RCP<PHX::DataLayout> >("Workset Scalar Data Layout")),
00063   J           (p.get<std::string>                   ("DetDefGrad Name"),
00064     p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00065   defgrad     (p.get<std::string>                   ("DefGrad Name"),
00066       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
00067     TResidual   (p.get<std::string>                   ("Residual Name"),
00068     p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
00069     haveSource  (p.get<bool>("Have Source")),
00070     haveConvection(false),
00071     haveAbsorption  (p.get<bool>("Have Absorption")),
00072     haverhoCp(false)
00073   {
00074     if (p.isType<bool>("Disable Transient"))
00075       enableTransient = !p.get<bool>("Disable Transient");
00076     else enableTransient = true;
00077 
00078     this->addDependentField(stabParameter);
00079     this->addDependentField(deltaTime);
00080     this->addDependentField(weights);
00081     this->addDependentField(wBF);
00082     this->addDependentField(porePressure);
00083     this->addDependentField(ThermalCond);
00084     this->addDependentField(kcPermeability);
00085     this->addDependentField(porosity);
00086     this->addDependentField(biotCoefficient);
00087     this->addDependentField(biotModulus);
00088     this->addDependentField(young_modulus_);
00089     this->addDependentField(poissons_ratio_);
00090     this->addDependentField(densityPoreFluid);
00091 
00092     this->addDependentField(Temp);
00093     this->addDependentField(RefTemp);
00094     this->addDependentField(TGrad);
00095     this->addDependentField(TempGrad);
00096     this->addDependentField(wGradBF);
00097 
00098     this->addDependentField(J);
00099     this->addDependentField(alphaMixture);
00100     this->addDependentField(alphaSkeleton);
00101     this->addDependentField(alphaPoreFluid);
00102     this->addDependentField(defgrad);
00103     this->addEvaluatedField(TResidual);
00104 
00105     if (haveSource) this->addDependentField(Source);
00106       if (haveAbsorption) {
00107         Absorption = PHX::MDField<ScalarT,Cell,QuadPoint>(
00108                 p.get<std::string>("Absorption Name"),
00109                 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"));
00110         this->addDependentField(Absorption);
00111       }
00112 /*
00113     Teuchos::RCP<PHX::DataLayout> vector_dl =
00114       p.get< Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout");
00115     std::vector<PHX::DataLayout::size_type> dims;
00116     vector_dl->dimensions(dims);
00117 */
00118     Teuchos::RCP<PHX::DataLayout> vector_dl =
00119       p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00120     std::vector<PHX::DataLayout::size_type> dims;
00121     vector_dl->dimensions(dims);
00122     numQPs  = dims[1];
00123     numDims = dims[2];
00124 
00125     Teuchos::RCP<PHX::DataLayout> node_dl =
00126       p.get< Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout");
00127     std::vector<PHX::DataLayout::size_type> ndims;
00128     node_dl->dimensions(ndims);
00129     worksetSize = dims[0];
00130     numNodes = dims[1];
00131 
00132     // Get data from previous converged time step
00133     porePressureName = p.get<std::string>("QP Pore Pressure Name")+"_old";
00134     JName =p.get<std::string>("DetDefGrad Name")+"_old";
00135     TempName =p.get<std::string>("QP Temperature Name")+"_old";
00136 
00137     // Works space FCs
00138     C.resize(worksetSize, numQPs, numDims, numDims);
00139     Cinv.resize(worksetSize, numQPs, numDims, numDims);
00140     F_inv.resize(worksetSize, numQPs, numDims, numDims);
00141     F_invT.resize(worksetSize, numQPs, numDims, numDims);
00142     JF_invT.resize(worksetSize, numQPs, numDims, numDims);
00143     KJF_invT.resize(worksetSize, numQPs, numDims, numDims);
00144     Kref.resize(worksetSize, numQPs, numDims, numDims);
00145 
00146     // Allocate workspace
00147     flux.resize(dims[0], numQPs, numDims);
00148     fgravity.resize(dims[0], numQPs, numDims);
00149     fluxdt.resize(dims[0], numQPs, numDims);
00150     pterm.resize(dims[0], numQPs);
00151     Tterm.resize(dims[0], numQPs);
00152 
00153     tpterm.resize(dims[0], numNodes, numQPs);
00154 
00155     if (haveAbsorption)  aterm.resize(dims[0], numQPs);
00156 
00157     convectionVels = Teuchos::getArrayFromStringParameter<double> (p,
00158                    "Convection Velocity", numDims, false);
00159     if (p.isType<std::string>("Convection Velocity")) {
00160       convectionVels = Teuchos::getArrayFromStringParameter<double> (p,
00161                      "Convection Velocity", numDims, false);
00162     }
00163     if (convectionVels.size()>0) {
00164       haveConvection = true;
00165       if (p.isType<bool>("Have Rho Cp"))
00166   haverhoCp = p.get<bool>("Have Rho Cp");
00167       if (haverhoCp) {
00168   PHX::MDField<ScalarT,Cell,QuadPoint> tmp(p.get<std::string>("Rho Cp Name"),
00169              p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"));
00170   rhoCp = tmp;
00171   this->addDependentField(rhoCp);
00172       }
00173     }
00174     this->setName("ThermoPoroPlasticityResidMass"+PHX::TypeString<EvalT>::value);
00175   }
00176 
00177   //**********************************************************************
00178   template<typename EvalT, typename Traits>
00179   void ThermoPoroPlasticityResidMass<EvalT, Traits>::
00180   postRegistrationSetup(typename Traits::SetupData d,
00181       PHX::FieldManager<Traits>& fm)
00182   {
00183   this->utils.setFieldData(stabParameter,fm);
00184     this->utils.setFieldData(Temp,fm);
00185     this->utils.setFieldData(RefTemp,fm);
00186   this->utils.setFieldData(deltaTime,fm);
00187   this->utils.setFieldData(weights,fm);
00188     this->utils.setFieldData(wBF,fm);
00189     this->utils.setFieldData(porePressure,fm);
00190     this->utils.setFieldData(ThermalCond,fm);
00191     this->utils.setFieldData(kcPermeability,fm);
00192     this->utils.setFieldData(porosity,fm);
00193     this->utils.setFieldData(alphaMixture,fm);
00194     this->utils.setFieldData(biotCoefficient,fm);
00195     this->utils.setFieldData(biotModulus,fm);
00196     this->utils.setFieldData(young_modulus_,fm);
00197     this->utils.setFieldData(poissons_ratio_,fm);
00198     this->utils.setFieldData(TGrad,fm);
00199     this->utils.setFieldData(TempGrad,fm);
00200     this->utils.setFieldData(wGradBF,fm);
00201     this->utils.setFieldData(J,fm);
00202     this->utils.setFieldData(defgrad,fm);
00203     this->utils.setFieldData(densityPoreFluid,fm);
00204     this->utils.setFieldData(alphaPoreFluid,fm);
00205     this->utils.setFieldData(alphaSkeleton,fm);
00206     this->utils.setFieldData(TResidual,fm);
00207     if (haveSource)  this->utils.setFieldData(Source,fm);
00208     if (haveAbsorption)  this->utils.setFieldData(Absorption,fm);
00209     if (haveConvection && haverhoCp)  this->utils.setFieldData(rhoCp,fm);
00210   }
00211 
00212 //**********************************************************************
00213 template<typename EvalT, typename Traits>
00214 void ThermoPoroPlasticityResidMass<EvalT, Traits>::
00215 evaluateFields(typename Traits::EvalData workset)
00216 {
00217   typedef Intrepid::FunctionSpaceTools FST;
00218   typedef Intrepid::RealSpaceTools<ScalarT> RST;
00219 
00220   Albany::MDArray porosityold = (*workset.stateArrayPtr)[porosityName];
00221   Albany::MDArray porePressureold = (*workset.stateArrayPtr)[porePressureName];
00222   Albany::MDArray Jold = (*workset.stateArrayPtr)[JName];
00223   Albany::MDArray Tempold = (*workset.stateArrayPtr)[TempName];
00224 
00225   ScalarT dTemperature(0.0);
00226   ScalarT dporePressure(0.0);
00227   ScalarT dJ(0.0);
00228 
00229   // Pore-Fluid Diffusion Term
00230 
00231    ScalarT dt = deltaTime(0);
00232 
00233    // Pull back permeability
00234    RST::inverse(F_inv, defgrad);
00235    RST::transpose(F_invT, F_inv);
00236    FST::scalarMultiplyDataData<ScalarT>(JF_invT, J, F_invT);
00237    FST::scalarMultiplyDataData<ScalarT>(KJF_invT, kcPermeability, JF_invT);
00238    FST::tensorMultiplyDataData<ScalarT>(Kref, F_inv, KJF_invT);
00239 
00240    /*
00241    // gravity or other potential term
00242      for (std::size_t cell=0; cell < workset.numCells; ++cell){
00243          for (std::size_t qp=0; qp < numQPs; ++qp) {
00244            for (std::size_t dim=0; dim <numDims; ++dim){
00245             fgravity(cell,qp, dim) = TGrad(cell,qp,dim);
00246          }
00247         fgravity(cell, qp, 1) -=  9.81*densityPoreFluid(cell, qp)*
00248                                       std::exp( porePressure(cell,qp)/biotModulus(cell,qp)-
00249                                            3.0* alphaPoreFluid(cell,qp)*
00250                                               (Temp(cell,qp) - RefTemp(cell,qp))); //assume g is 8.81
00251      }
00252    }
00253    */
00254 
00255    // Pore pressure gradient contribution
00256   FST::tensorMultiplyDataData<ScalarT> (flux, Kref, TGrad); // flux_i = k I_ij p_j
00257   // FST::tensorMultiplyDataData<ScalarT> (flux, Kref, fgravity); // flux_i = k I_ij p_j
00258 
00259    for (std::size_t cell=0; cell < workset.numCells; ++cell){
00260       for (std::size_t qp=0; qp < numQPs; ++qp) {
00261         for (std::size_t dim=0; dim < numDims; ++dim){
00262           fluxdt(cell, qp, dim) = -flux(cell,qp,dim)*dt;
00263         }
00264       }
00265   }
00266   FST::integrate<ScalarT>(TResidual, fluxdt, wGradBF, Intrepid::COMP_CPP, false); // "false" overwrites
00267 
00268   // Pore-fluid diffusion coupling.
00269   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00270     for (std::size_t node=0; node < numNodes; ++node) {
00271   //    TResidual(cell,node)=0.0;
00272       for (std::size_t qp=0; qp < numQPs; ++qp) {
00273 
00274             dJ = std::log(J(cell,qp)/Jold(cell,qp));
00275             dTemperature = Temp(cell,qp)-Tempold(cell,qp);
00276             dporePressure = porePressure(cell,qp)-porePressureold(cell, qp);
00277 
00278           // Volumetric Constraint Term
00279           TResidual(cell,node) -=  (biotCoefficient(cell, qp)*dJ
00280                                         - 3.0*alphaSkeleton(cell,qp)*J(cell,qp)*
00281                                         (Temp(cell,qp) - RefTemp(cell,qp))*dJ  )
00282                                     *wBF(cell, node, qp) ;
00283 
00284           // Pore-fluid Resistance Term
00285           TResidual(cell,node) -=  dporePressure/biotModulus(cell, qp)*
00286                                        wBF(cell, node, qp);
00287 
00288          // Thermal Expansion
00289          TResidual(cell,node) +=  3.0*dTemperature*alphaMixture(cell, qp)*
00290                                 wBF(cell, node, qp);
00291 
00292         }
00293       }
00294   }
00295 
00296   // Projection-based Stabilization
00297 
00298 // Penalty Term
00299   for (std::size_t cell=0; cell < workset.numCells; ++cell){
00300 
00301     porePbar = 0.0;
00302     Tempbar = 0.0;
00303     vol = 0.0;
00304     for (std::size_t qp=0; qp < numQPs; ++qp) {
00305 
00306     porePbar += weights(cell,qp)*(porePressure(cell,qp)-porePressureold(cell, qp));
00307     Tempbar += weights(cell,qp)*(Temp(cell,qp)-Tempold(cell,qp));
00308 
00309      vol  += weights(cell,qp);
00310      }
00311      porePbar /= vol;
00312      Tempbar /= vol;
00313 
00314      for (std::size_t qp=0; qp < numQPs; ++qp) {
00315        pterm(cell,qp) = porePbar;
00316        Tterm(cell,qp) = Tempbar;
00317    }
00318   }
00319   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00320     for (std::size_t node=0; node < numNodes; ++node) {
00321       for (std::size_t qp=0; qp < numQPs; ++qp) {
00322 
00323         dTemperature = Temp(cell,qp)-Tempold(cell,qp);
00324         dporePressure = porePressure(cell,qp)-porePressureold(cell, qp);
00325 
00326         shearModulus = young_modulus_(cell,qp)*0.5/(1.0+ poissons_ratio_(cell,qp));
00327         bulkModulus =  young_modulus_(cell,qp)/3.0/(1.0- 2.0*poissons_ratio_(cell,qp));
00328 
00329         safeFactor =(bulkModulus + 4.0/3.0*shearModulus +
00330                  biotCoefficient(cell,qp)*biotCoefficient(cell,qp)*biotModulus(cell,qp))/
00331                  (biotModulus(cell,qp)*(bulkModulus + 4.0/3.0*shearModulus));
00332 
00333         TResidual(cell,node) += (pterm(cell,qp)-dporePressure)
00334                                               *stabParameter(cell, qp)
00335                                               *safeFactor
00336                                                 *wBF(cell, node, qp);
00337 
00338               safeFactor = alphaMixture(cell,qp)*(bulkModulus + 4.0/3.0*shearModulus +
00339               biotCoefficient(cell,qp)*biotCoefficient(cell,qp)*biotModulus(cell,qp))/
00340               (bulkModulus + 4.0/3.0*shearModulus);
00341 
00342         TResidual(cell,node) += 3.0*(dTemperature - Tterm(cell,qp))
00343                                  *stabParameter(cell, qp)
00344                                  *safeFactor
00345                                  *wBF(cell, node, qp);
00346       }
00347     }
00348   }
00349 
00350 }
00351 //**********************************************************************
00352 }
00353 
00354
ThermoPoroPlasticityResidMass_Def.hpp
